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PREFACE 


This  is  the  first  volume  of  technical  reports  produced  for  the  project  ‘Laser 
Peening  for  Reliable  Fatigue  Life,’  Contract  No.  FA8650-04-D-3446-25.  The 
scope  of  this  volume  includes  the  simulation  of  the  Laser  Peening  process  in 
3D,  the  identification  of  the  key  process  parameters,  the  optimization  of  the 
residual  stress  profile  using  the  identified  key  parameters,  an  outline  of  the 
optimization  strategies  employed,  and  the  modifications  made  to  Particle 
Swarm  Optimization  to  yield  multiple  optima. 

Additional  volumes  cover  various  other  aspects  of  this  project.  The  next  of 
these  volumes  covers  the  experimental  validation  of  the  fatigue  life  model, 
which  empirically  demonstrates  the  accuracy  of  the  fatigue  life  model  used  in 
this  work.  The  error  between  the  model  and  reality  is  addressed  and  the 
sources  of  this  error  are  investigated  and  amended  where  possible.  These 
sources  include  material  parameters  used  in  the  Finite  Element  Modeling, 
accuracy  of  strain-based  fatigue  life  model  parameters,  and  residual  stress 
relaxation  effects.  The  third  topic  covers  the  inverse  optimization  of  high 
strain-rate  material  properties  from  experimental  material  data.  The 
constructed  model  is  then  validated  and  quantification  of  uncertainty 
developed.  The  fourth  topic  addresses  miscellaneous  effects  of  Laser  Peened 
residual  stress  fields  on  fatigue  crack  growth  and  machined  edges.  In 
addition.  Laser  Peening  residual  stress  fields  are  adapted  from  fiat  to  curved 
plates.  Lastly,  the  development  of  model-form  and  parametric  uncertainty  are 
developed  as  well  as  the  inclusion  of  both  types  simultaneously.  This 
uncertainty  quantification  is  developed  from  the  material  model  analysis 
conducted  and  incorporates  results  from  the  experimental  work  and  realistic 
variation  within  Laser  Peening  input  parameters. 
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SUMMARY 


Laser  peening  (LP)  is  a  surface  enhancement  technique  that  has  been  applied 
to  improve  fatigue  and  corrosion  properties  of  metals.  The  ability  to  use  a  high 
energy  laser  pulse  to  generate  shock  waves,  inducing  a  compressive  residual 
stress  field  in  metallic  materials,  has  applications  in  multiple  fields  such  as 
turbomachinery,  airframe  structures,  and  medical  appliances.  In  the  past,  re¬ 
searchers  have  investigated  the  effects  of  LP  parameters  experimentally  and 
performed  a  limited  number  of  simulations  on  simple  geometries.  However, 
monitoring  the  dynamic,  intricate  relationships  of  peened  materials  experi¬ 
mentally  is  time  consuming,  expensive,  and  challenging. 

With  increasing  applications  of  LP  on  complex  geometries,  these  limited  ex¬ 
perimental  and  simulation  capabilities  are  not  sufficient  for  an  effective  LP 
process  design.  Due  to  high  speed,  dynamic  process  parameters,  it  is  diffi¬ 
cult  to  achieve  a  consistent  residual  stress  field  in  each  treatment  and  constrain 
detrimental  effects.  With  increased  computer  speed  as  well  as  increased  so¬ 
phistication  in  non-linear  finite  element  analysis  software,  it  is  now  possible  to 
develop  simulations  that  can  consider  several  LP  parameters. 

In  this  research,  a  finite  element  simulation  capability  of  the  LP  process  is 
developed.  These  simulations  are  validated  with  the  available  experimental 
results.  Based  on  the  validated  model,  simplifications  to  complex  models  are 
developed.  These  models  include  quarter  symmetric  3D  model,  a  cylindrical 
coupon,  a  parametric  plate,  and  a  bending  coupon  model.  The  developed  mod¬ 
els  can  perform  simulations  incorporating  the  LP  process  parameters,  such 
as  pressure  pulse  properties,  spot  properties,  number  of  shots,  locations,  se¬ 
quences,  overlapping  configurations,  and  complex  geometries.  These  models 
are  employed  in  parametric  investigations  and  residual  stress  profile  optimiza¬ 
tion  at  single  and  multiple  locations. 

In  parametric  investigations,  quarter  symmetric  3D  model  is  used  to  investigate 
temporal  variations  of  pressure  pulse,  pressure  magnitude,  and  shot  shape  and 
size.  The  LP  optimization  problem  is  divided  into  two  parts:  single  and  mul¬ 
tiple  locations  peening  optimization.  The  single-location  peening  optimiza¬ 
tion  problems  have  mixed  design  variables  and  multiple  optimal  solutions.  In 
the  optimization  literature,  many  researchers  have  solved  problems  involving 
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mixed  variables  or  multiple  optima,  but  it  is  difficult  to  find  multiple  solutions 
for  mixed- variable  problems.  A  mixed- variable  Niche  Particle  Swarm  Opti¬ 
mization  (MNPSO)  is  proposed  that  incorporates  a  mixed-variable  handling 
technique  and  a  niching  technique  to  solve  the  problem. 

Designing  an  optimal  residual  stress  profile  for  multiple-location  peening  is  a 
challenging  task  due  to  the  computational  cost  and  the  nonlinear  behavior  of 
LP.  A  Progressive  Multifidelity  Optimization  Strategy  (PMOS)  is  proposed  to 
solve  the  problem.  The  three-stage  PMOS,  combines  low-  and  high-  fidelity 
simulations  and  respective  surrogate  models  and  a  mixed-variable  handling 
strategy.  This  strategy  employs  comparatively  low  computational-intensity 
models  in  the  first  two  stages  to  locate  the  design  space  that  may  contain  the 
optimal  solution.  The  third  stage  employs  high  fidelity  simulation  and  surro¬ 
gate  models  to  determine  the  optimal  solution.  The  overall  objective  of  this 
research  is  to  employ  finite  element  simulations  and  effective  optimization 
techniques  to  achieve  optimal  residual  stress  fields. 
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1.0  INTRODUCTION 


Material  is  subjected  to  degradation  from  the  moment  it  is  prepared.  Mechan¬ 
ical,  chemical,  or  thermal  loading  continues  to  degrade  the  material  through¬ 
out  its  lifetime.  Damage  by  fretting,  fatigue,  corrosion,  and  wear  are  a  few 
examples  of  the  degradation  that  can  inflict  a  cost  penalty  on  a  material  by 
requiring  design  changes  to  accommodate  degradation.  Material  damage  can 
also  cause  loss  of  performance.  Most  types  of  material  degradation  tend  to 
originate  at  the  surface  of  a  component.  The  objective  of  surface  enhancement 
techniques  is  to  inhibit  degradation  by  altering  material  properties  at  the  sur¬ 
face  and  to  maintain  the  required  level  of  performance  without  imposing  high 
cost  implications.  However,  it  is  an  intriguing,  but  challenging  to  achieve  the 
objective  by  tailoring  the  surface  enhancement  techniques.  In  this  research, 
a  simulation-based  strategy  is  employed  on  a  selected  surface  enhancement 
technique  to  engineer  optimal  surface  properties. 

The  main  focus  of  this  research  is  the  surface  enhancement  technique  called 
Laser  Peening  (LP).  In  the  LP  process,  favorable  residual  stresses  are  induced 
on  a  surface  to  improve  fatigue  and  fretting  properties  of  metals.  In  the  litera¬ 
ture,  experimental  [1-71  and  simulation  [8-111  work  have  been  performed  to 
drive  maximum  benefits  from  the  process.  However,  time  consuming  and  ex¬ 
pensive  experiments  and  limited  simulations  on  simple  geometries  are  not  suf¬ 
ficient  for  optimal  LP  process  design.  A  comprehensive  procedure  is  required 
that  can  not  only  perform  simulations  but  also  can  be  employed  in  obtaining 
optimal  process  parameters. 

The  overall  goal  of  this  section  is  to  provide  a  brief  introduction  of  surface 
enhancement  techniques  and  to  present  a  detailed  description  of  a  specific  pro¬ 
cess  called  Laser  Peening.  This  section  begins  with  Section  1.1,  an  introduc¬ 
tion  to  surface  enhancement  techniques  such  as  shot  peening,  low  plasticity 
burnishing,  and  waterjet  peening,  and  their  use  in  manufacturing  processes. 
Section  1.2  presents  an  introduction  to  the  LP  process.  The  next  three  sections 
present  the  laser  generation  mechanism,  component  preparation,  and  shock 
wave  propagation.  In  Section  1.3  includes,  a  preliminary  explanation  of  a 
laser  generation  mechanism.  Section  1.4  presents  component  surface  prepara¬ 
tion  details.  After  the  component  is  prepared,  a  laser  is  fired  on  the  surface. 
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Section  1.5  deals  with  the  shock  wave  propagation  and  residual  stress  genera¬ 
tion  in  a  component.  Section  1.6  discusses  the  advantages  and  disadvantages 
of  LP  compared  to  other  mechanical  surface  enhancement  techniques.  Section 
1.7  presents  the  organization  of  the  document,  and  the  final  section  summa¬ 
rizes  the  section. 

1.1  Surface  Enhancement  Processes 

Connecting  rods  in  automobiles,  turbine  blades  in  aircraft  engines,  and  lugs  in 
aircraft  structures  are  a  few  components  that  experience  material  damage  due 
to  fretting  and  fatigue  loading  conditions.  These  type  of  loads  cause  failure  by 
initializing  and  propagating  cracks.  In  the  case  of  turbine  engines,  foreign  ob¬ 
ject  damages  result  in  early  crack  initiation  and  growth,  causing  a  decrease  in 
fatigue  performance.  The  possibility  of  such  failures  leads  to  increased  inspec¬ 
tion  and  maintenance  expenses  because  the  failure  of  one  critical  component 
limits  the  operational  time  and  service  life  of  the  entire  system.  Fatigue  can 
cause  catastrophic  failure,  leading  to  structural  and  automobile  failure,  aircraft 
losses,  and,  at  worst,  loss  of  life.  Apart  from  standard  design  requirements 
and  inspections,  many  material  treatments  have  been  developed  to  mitigate 
such  failures.  Surface  enhancement  techniques  treat  the  intended  surfaces  to 
improve  the  desired  surface  properties.  These  techniques  are  becoming  an  in¬ 
tegral  part  of  manufacturing  processes  because  the  behavior  of  manufactured 
parts  is  dependent  on  the  structure  and  properties  both  of  bulk  material  and  of 
the  surface. 

Every  metallic  material  consists  of  a  micro- structure  of  small  crystals  called 
grains  or  crystallites.  The  properties  of  these  grains  (i.e.  grain  size,  orienta¬ 
tion,  and  composition)  determine  the  overall  behavior  of  the  material.  Mate¬ 
rial  treatment  (and/or  surface  enhancement  techniques)  provides  an  efficient 
way  to  influence  the  properties  of  the  metal  by  controlled  processes.  These 
material  treatments  are  generally  classified  into  mechanical,  thermal,  thermo¬ 
mechanical,  and  thermo-chemical  treatments.  In  the  case  of  mechanical  treat¬ 
ments,  some  form  of  mechanical  force  is  employed  to  modify  the  properties 
of  the  material  to  the  advantage  of  an  application.  Common  mechanical  treat¬ 
ment  techniques  are  shot  peening,  low  plasticity  burnishing,  water  jet  peening, 
cold  rolling,  and  case  hardening. 
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In  the  case  of  thermal  treatments,  some  form  of  temperature  controls  are  used 
to  modify  material  properties.  Common  thermal  treatment  techniques  are 
hardening,  tempering,  and  annealing.  Hot  rolling  and  plating  are  examples 
of  thermo-mechanical  treatment,  a  category  that  combines  the  effect  of  both 
mechanical  and  thermal  treatments.  In  the  case  of  thermo-chemical  treatments, 
a  metallic  or  non-metallic  coating  is  applied  on  the  surface.  Painting,  oxidiz¬ 
ing,  and  vapor  deposition  are  a  few  examples  of  thermo-chemical  treatments. 
Depending  upon  the  usage,  a  few  of  these  processes  are  also  referred  to  as 
surface  enhancement  techniques.  Surface  enhancement  techniques  have  sig¬ 
nificant  industrial  applications.  The  following  section  talks  about  mechanical 
surface  enhancement  techniques. 

1.1.1  Shot  Peening 

Shot  peening  (SP)  is  a  traditional  surface  enhancement  technique  that  has  been 
used  for  more  than  six  decades.  SP  was  discovered  and  patented  by  Lockheed 
Aircraft  Corporation  in  Burbank,  California  in  the  late  1940’s.  For  example, 
Boeing  Commercial  Aircraft  Company  have  standard  setup  to  peening  aircraft 
components  [121.  Many  other  companies  are  using  SP,  or  a  modified  form. 
Glass  Bead  Peening  (GBP),  to  improve  fatigue  life  of  mechanical  components. 
In  the  SP  process,  small,  spherical  metal  or  ceramic  balls  are  bombarded  on 
the  surface  to  be  peened.  In  the  GBP  process,  glass  balls  are  used  instead  of 
metallic  balls.  This  bombardment  creates  elastic  and  plastic  deformation.  This 
deformation  tends  to  produce  compressive  stress  on  the  surface  and  tensile 
stress  in  the  interior.  The  schematic  of  the  SP  process  is  shown  in  Figure 
1.  SP  results  in  a  small  depth  of  induced  compressive  residual  stress.  This 
process  results  in  a  roughened  surface  due  to  the  physical  contact  between  the 
bombarding  ball  and  the  peened  component.  This  effect  is  more  severe  in  soft 
materials,  and  may  not  be  an  ideal  process  if  the  surface  finish  is  an  important 
factor  for  a  product.  With  proper  parameters,  the  SP  technique  [13-171  is 
used  to  induce  compressive  residual  stresses  in  components.  The  physical 
principle  of  changing  the  mechanical  properties  is  approximately  the  same  for 
most  mechanical  surface  enhancement  techniques.  The  principle  is  that  elastic 
and  plastic  deformation  result  in  residual  stress  generation.  This  generation  is 
studied  in  detail  in  Section  1.2. 
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Figure  1 :  Schematic  of  shot  peening  process 


1.1.2  Low  Plasticity  Burnishing 

Low  Plasticity  Burnishing  (LPB)  is  a  recently  developed  (1998)  surface  en¬ 
hancement  technique.  The  LPB  process  [18, 191  consists  of  rolling  a  high- 
modulus  ball  or  roller  over  the  intended  surface  of  the  component  (Figure  2). 
Depending  upon  the  controlling  parameters,  this  process  alters  the  mechanical 
behavior  of  the  material  by  low-cold  working,  just  like  SR  This  is  a  low-cost, 
easy-to  set-up  process  that  provides  a  better  surface  finish  compared  to  SR  An 
advantage  of  this  process  is  that  it  can  be  carried  out  in  any  numerically  con¬ 
trolled  apparatus,  including  a  lathe,  mill,  or  CNC  machine.  The  surface  finish 
obtained  depends  upon  the  finish  of  the  ball  used  in  LPB  process;  the  ball  is  a 
wear-prone  component.  Just  like  SP,  this  process  produces  elastic  and  plastic 
deformation  by  contact  between  the  burnishing  ball  and  the  surface.  A  disad¬ 
vantage  of  the  LPB  process  is  that  it  is  difficult  to  apply  to  curved  or  complex 
geometries. 
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Normal  force 


1.1.3  Waterjet  Peening 

High  pressure  waterjets  have  been  studied  to  understand  the  mechanisms  as¬ 
sociated  with  the  jet  material  interface  and  to  apply  it  to  multiple  processes 
such  as  cleaning,  cutting,  and  paint  removal.  The  application  of  an  ultra-high- 
pressure  waterjet  as  a  surface  enhancement  technique  is  relatively  new  (2000) 
compared  to  SP  and  LPB.  Waterjet  Peening  or  Water  Peening  (WP)  [20, 21] 
is  similar  to  SP  except  that  it  uses  high  pressure  droplets  that  disintegrate  in 
the  waterjet  flow  fleld  instead  of  solid  shots  [22].  Compared  to  other  avaflable 
techniques,  the  tower  cost,  the  absence  of  heat  that  coutd  affect  the  region,  and 
a  ctean  surface  are  the  major  advantages  of  the  process.  The  WP  process  is  a 
physicatly  complex  technique  that  is  difficult  to  model  and  requires  extensive 
further  research. 

1.2  Laser  Peening 

High-power  Light  Ampliflcation  by  Simulated  Emission  of  Radiation  (LASER) 
is  used  in  many  manufacturing  processes  of  aircraft  engine  and  airframe  struc¬ 
tures,  such  as  laser  drilling,  welding,  and  glazing.  LP  is  a  surface  enhancement 
technique  for  metallic  components  that  uses  the  shock  wave  generation  capa¬ 
bility  of  lasers  to  induce  favorable  residual  stresses.  Apart  from  improving 
fatigue  life,  LP  has  also  been  applied  to  improve  fretting  and  corrosion  prop¬ 
erties. 
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The  ability  of  a  pulsed  laser  beam  to  generate  shock  waves  was  first  recog¬ 
nized  and  explored  in  the  early  1960’s.  Initial  facilities  were  developed  and 
feasibility  studies  were  performed  at  Battelle  Laboratories  in  Columbus,  OH. 
Also  at  Battelle,  researchers  succeeded  in  the  application  of  LP  to  enhance 
fatigue  properties  of  fastened  joints  for  aeronautical  applications.  After  this 
success,  further  research  was  performed  to  examine  various  fundamental  prin¬ 
ciples  such  as  the  confined  interaction  mode,  the  dielectric  breakdown  factor, 
and  analytical  modeling  of  the  physical  process  involved.  As  a  result  of  these 
numerous  efforts,  LP  is  emerging  as  an  alternative  and  complementary  process 
to  conventional  peening  processes. 

Compared  with  traditional  surface  enhancement  techniques  such  as  SP  and 
LPB,  a  higher  depth  of  compressive  residual  stress  and  a  lower  cold- work  am¬ 
plitude  are  obtained  from  the  LP  treatment.  Among  the  major  advantages  of 
the  LP  process,  four  important  ones  are  the  absence  of  contact  (prevents  se¬ 
vere  relief  of  residual  stresses  during  cyclic  loading),  better  surface  finish  than 
shot  peening,  precise  parameter  control,  and  the  ability  to  peen  complicated 
geometries. 

1.2.1  Laser  Generation 

To  generate  a  laser  [23,241,  the  atoms  and  molecules  of  a  crystal,  gas,  or  liquid 
are  excited  so  that  more  of  these  are  at  higher  energy  levels  than  are  at  lower 
energy  levels.  When  a  photon  of  the  frequency  that  corresponds  to  the  energy 
difference  between  the  excited  and  lower  states  strikes  an  excited  atom,  this 
process  causes  atoms  to  fall  back  to  a  lower  energy  state  and  to  emit  a  second 
photon  of  the  same  or  a  proportional  frequency,  in  phase  with  and  in  the  same 
direction  as  the  striking  photon.  This  process  is  called  stimulated  emission,  in 
which  an  atom  is  stimulated  to  produce  a  second  photon.  The  striking  and  the 
emitted  photon  may  stimulate  further  atoms  to  emit  photons,  all  with  the  same 
frequency  and  phase.  This  process  is  a  rapid  chain  reaction  that  produces  a 
sudden  burst  of  coherent  radiation  as  all  the  atoms  discharge. 

Based  on  different  laser  parameters,  there  are  many  varieties  of  lasers  avail¬ 
able.  The  intensity,  duration  (mid-span),  and  wave  length  of  the  laser  are 
significant  properties.  In  the  case  of  the  LP  process,  these  properties  should 


be  such  that  the  laser  produces  the  required  volume,  depth,  and  magnitude 
of  residual  stress  in  the  peened  material.  The  LP  process  requires  a  laser 
with  high  power  (100  to  1000  watts),  moderate  frequency  (10  to  50  shots  per 
minute),  and  high  laser  pulse  energy  (100  Jules)  for  suitable  industrial  applica¬ 
tions.  A  very  high  repetition  rate  may  not  be  needed  because  in  the  LP  process, 
repetition  is  limited  setup  for  each  shot.  Consequently,  the  laser  peening  speed 
is  limited  by  component  preparation  time  rather  than  laser  frequency. 

1.2.2  Component  Preparation 

The  ability  of  a  laser  to  create  stresses  in  a  metal  component  was  recognized 
far  in  advance  of  its  practical  application.  In  the  process  of  bringing  LP  to  the 
level  of  industrial  application,  the  use  of  opaque  and  transparent  overlays  have 
played  a  significant  role.  To  generate  plasma,  the  component  surface  must 
be  painted  with  an  opaque  overlay.  The  sudden  expansion  of  plasma  causes 
a  pressure  pulse  in  the  metal.  The  pressure  pulse  peak  and  mid-span  can  be 
increased  if  the  surface  is  confined  with  a  transparent  overlay.  The  role  of 
overlays  is  further  discussed  in  the  following  sections. 

Opaque  Overlay 

When  the  laser  beam  is  directed  onto  the  surface,  it  passes  through  the  trans¬ 
parent  overlay  and  strikes  the  opaque  overlay.  This  produces  a  high  tempera¬ 
ture  of  the  order  of  10,000®F  during  the  plasma  creation  stage.  Direct  contact 
of  a  metal  surface  with  the  high  temperature  plasma  will,  in  most  cases,  form 
a  thin  melted  layer.  Therefore  at  the  nascent  stage  of  the  process,  the  metal 
component  is  covered  with  a  protective  coating  to  avoid  the  thermal  effects  of 
such  high  temperatures. 

However,  the  role  of  this  coating  has  changed  over  time.  At  present  time, 
the  opaque  overlay  has  two  purposes:  augmenting  the  creation  of  plasma  and 
providing  protection  of  the  component  surface  to  be  peened.  For  the  first  pur¬ 
pose,  the  opaque  overlay  supplies  material  for  plasma  creation.  For  the  sec¬ 
ond  purpose,  it  provides  protection  to  the  component  surface  from  high  tem¬ 
perature.  Because  of  this  protection,  the  material  undergoes  reduced  thermal 
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and  micro- structural  changes.  If  allowed,  these  thermal  and  micro- structural 
changes  can  alleviate  the  required  effects  of  LP.  The  important  factors  in  se¬ 
lecting  an  opaque  overlay  material  for  a  particular  material  surface  are  the 
cost,  ease  of  application  and  removal,  layer  thickness,  interaction  with  laser, 
opaque  integrity  for  multiple  shots,  applicability  over  complex  geometries,  dis¬ 
posal  of  used  opaque  material,  and  repeatability.  Among  all  opaque  coatings 
black  paint,  aluminum,  zinc  or  copper  and  organic  coating  are  used  commonly. 
These  are  preferred  because  of  the  ease  of  handling  and  disposal. 


Transparent  Overlay 

When  a  dielectric  material,  transparent  to  laser  light,  is  placed  over  the  opaque 
overlay,  plasma  created  by  LP  is  trapped  between  the  component  and  the  trans¬ 
parent  overlay.  This  trapping  increases  the  pressure  pulse  magnitude  by  factors 
of  5  (for  the  peak  pressure  level)  and  3  (for  the  pulse  mid- span/duration)  com¬ 
pared  to  a  direct  ablation  mode. 

In  industrial  applications  where  a  large  number  of  shots  is  needed  in  less  time, 
water  is  used  most  often  for  the  overlay  because  of  ease  of  availability  and 
handling.  However,  solid  overlays  of  glass,  fused  quartz,  and  acrylic  have  also 
been  used.  As  with  the  opaque  overlay,  the  cost,  speed  of  control,  ease  of 
application  and  removal,  and  ease  of  clean-up  and  disposal  are  major  factors 
in  selecting  the  transparent  overlay. 

Along  with  these  factors,  dielectric  breakdown  thresholds,  degree  of  trans¬ 
parency,  acoustic  impedance,  and  control  of  thickness  also  affect  selection  of 
the  overlay.  Proper  selection  of  the  transparent  overlay  is  important  because 
it  ultimately  affects  the  peak  and  duration  of  the  shock  wave  in  the  material. 
It  is  important  to  remember  that  in  the  industrial  application  of  opaque  and 
transparent  overlays,  a  traditional  manual  application  of  overlays  may  not  be 
most  suitable.  Instead,  a  liquid  jet  is  used  for  efficient  processing,  as  shown  in 
Figure  3.  These  jets  and  the  laser  beam  must  be  properly  sequenced  to  produce 
the  desired  results.  Proper  positioning  of  the  component  is  also  necessary  for 
successful  completion  of  the  process. 
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Figure  3:  Schematic  of  Laser  Peening  process 
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1.2.3  Residual  Stresses 


The  component  surface  to  be  peened  is  coated  with  a  opaque  material  and  cov¬ 
ered  with  a  transparent  overlay.  When  the  component  is  ready  for  peening,  a 
high  energy  laser  is  fired  on  the  prepared  surface.  The  laser  passes  through  the 
transparent  layer  and  causes  the  opaque  overlay  to  vaporise  and  then  convert 
into  plasma  when  it  absorbs  more  energy.  As  the  vapor  continues  to  absorb  the 
laser  energy,  it  expands.  As  the  confined  plasma  expands,  it  generates  a  high 
pressure,  short  duration  shock  wave  that  travels  into  the  base  material  and  the 
water.  The  presence  of  water  tends  to  confine  the  energy  and  increase  the  in¬ 
tensity  of  the  pressure  pulse  in  the  component.  Depending  upon  the  magnitude 
of  the  pressure  wave  in  the  base  metal,  this  can  cause  elastic  and  plastic  defor¬ 
mation.  This  deformation  generates  compressive  and  tensile  residual  stresses. 
The  mechanics  of  the  residual  stress  generation  process  are  explained  below. 

Shock  Wave  Propagation 

Depending  upon  its  magnitude,  an  LP  pressure  pulse  can  cause  plastic  de¬ 
formation  in  the  top  region  of  a  component.  As  the  shock  wave  progresses 
through  the  component,  its  magnitude  is  reduced  according  to  the  attenuation 
rate.  After  a  certain  depth,  the  shock  wave  magnitude  is  below  the  proportional 
limit  and  can  only  cause  elastic  deformation.  The  final  outcome  of  the  process 
is  a  residual  stress  field  in  the  component.  The  top  region  of  the  component 
tends  to  have  compressive  stresses,  followed  by  a  tensile  region  beneath  it. 

In  general,  the  study  of  stress  propagation  in  a  non-linear  medium  is  extremely 
difficult.  This  section  describes  the  physical  behavior  of  the  material  under  LP 
loading.  Material  behaves  differently  under  different  strain-rates.  In  the  case 
of  LP  loading,  a  typical  strain-rate  is  up  to  10^/^',  which  is  extremely  higher 
than  conventional  strain-rates  of  10“^/^'.  This  higher  strain-rate  is  referred  to 
as  shock  wave  or  impact  loading.  Material  behaves  differently  under  static 
(or  low  strain-rate)  and  dynamic  (high  strain-rate  or  shock  wave,  >  10^/.?) 
loadings  because  higher  stress  is  necessary  to  reach  the  yield  limit  in  an  impact 
loading  than  in  case  of  slow  loading. 

The  dynamic  yield  strength  of  a  material  is  a  function  of  strain,  strain-rate. 
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Figure  4:  Shock  wave  in  a  material 


and  temperature.  When  a  shock  wave  of  a  magnitude  above  the  dynamic  yield 
strength  limit  travels  through  a  material,  it  takes  a  certain  shape  (Figure  4) 
based  on  the  pressure  magnitude  and  material  properties  [25,261.  An  important 
parameter  for  this  shape  is  the  Hugoniot  Elastic  Limit  (HEL),  which  is  defined 
as  the  axial  stress  required  for  plastic  deformation  in  a  uniaxial  strain  state. 
The  relationship  between  the  HEL  (Ghel)  and  the  dynamic  yield  stress  (Oyd) 
is 


/l-2v\ 
Oyd  =  GHEL  1  ^  ) 


(1) 


where  V  is  Poission’s  ratio. 

LP  pressure,  when  applied,  generates  two  stress  waves,  elastic  and  plastic,  in 
the  material  [25,261.  According  to  stress  wave  propagation  theory,  the  speed 
of  the  elastic  (Vg)  and  plastic  (Vp)  waves  are  given  in  Eqs.  2  and  3,  respectively. 
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1/2 


(2) 


K  = 


£(l-v 


L(1+v)(1-2v)pJ 


L3(l-2v)p 


1/2 


(3) 


where  p  denotes  mass  density  and  E  denotes  Young’s  modulus.  The  speed 
of  the  elastic  wave  (Eq.  2)  is  faster  than  the  plastic  wave  (Eq.  3).  When  the 
applied  pressure  is  removed,  an  unloading  (release)  wave  travels  in  the  same 
direction.  If  the  loading  is  compressive,  then  the  release  wave  is  always  tensile. 

The  LP-generated  elastic  wave  travels  until  it  reaches  the  boundary  of  the  ma¬ 
terial  and  then  reflects  back.  The  speed  of  the  release  wave  is  greater  than  a 
plastic  wave,  so  depending  upon  the  loading  time,  it  is  possible  for  the  waves  to 
meet.  The  release  wave  reduces  the  amplitude  of  the  plastic  wave  and  reflects 
back  toward  the  starting  surface.  In  the  meantime,  the  first  elastic  wave,  which 
reflected  back  from  the  boundary  of  the  component,  also  meets  the  plastic 
wave.  There  is  a  complicated  wave  interaction  that  takes  place  between  these 
waves.  The  timing  of  the  interactions  between  these  waves  can  be  different  de¬ 
pending  upon  many  factors  such  as  the  material  properties,  the  relative  speed 
of  the  elastic  and  plastic  waves,  and  the  total  length  of  the  component.  This 
results  in  a  complex  distribution  of  elastic  and  plastic  strains  and  stresses,  en¬ 
suring  complex  dynamic  response  in  the  laser-peened  material.  Therefore,  the 
dynamic  responses  of  the  laser  peened  material  are  complex.  To  understand 
the  above  phenomenon  and  predict  the  final  residual  stress  profile,  an  effective 
simulation  methodology  is  required. 


Residual  Stress  Generation 

The  complicated  interactions  between  different  waves  inside  the  material  pro¬ 
duce  compressive  and  tensile  stress  regions.  The  compressive  stress  dominates 
in  the  top  region  of  the  component.  Below  this  compressive  stress  region,  there 
is  a  tensile  stress  region.  The  magnitude  of  the  tensile  stress  is  less  than  that 
of  the  compressive  stress.  However,  the  volume  of  the  tensile  stress  region  is 
higher  than  the  volume  of  the  compressive  region.  Depending  upon  compo¬ 
nent  geometry,  material  properties,  and  pressure  pulse  properties,  the  tensile 
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region  may  be  in  between  the  two  compressive  regions.  Depending  upon  the 
material  properties  and  the  component  dimensions,  the  compressive  region  at 
the  bottom  may  or  may  not  be  present.  The  presence  of  compressive  stress 
tends  to  reduce  the  chances  of  crack  initiation  and  growth.  The  tensile  stress 
region  is  generated  to  balance  the  compressive  stress  for  equilibrium  in  the 
absence  of  external  force. 

The  physics  behind  this  complex  formation  of  stress  is  explained  below.  Figure 
5  shows  a  schematic  of  the  residual  stress  generation  and  distribution.  Above 
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Figure  5:  Schematic  of  residual  stress  generation 

the  elastic  limit,  pressure  causes  the  surface  layer  of  the  target  to  expand  be¬ 
yond  the  elastic  limit.  Depending  upon  the  magnitude  of  the  pressure,  the 
layer  below  the  surface  layer  also  expands,  but  within  the  elastic  limit.  This 
expansion  is  in  the  normal  direction  to  the  applied  pressure.  According  to  the 
attenuation  rate,  the  expansion  reduces  along  the  depth  in  the  target.  After  a 
certain  depth,  the  expansion  is  elastic. 

The  elastically  deformed  part  of  the  component  attempts  to  return  back  to  its 
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original  shape,  but  plastic  deformation  is  permanent.  Therefore,  as  shown  in 
Figure  5,  the  top  layer  of  the  target  is  subjected  to  compressive  stress,  and  the 
area  below  this  layer  is  subjected  to  tensile  stress,  maintaining  equilibrium  in 
the  target  without  external  force.  An  increase  in  the  compressive  stress  volume 
or  magnitude  will  always  come  with  an  increase  in  the  tensile  stress  volume 
or  magnitude.  An  advantage  of  the  LP  process  is  a  higher  depth  of  compres¬ 
sive  stress  as  compare  to  other  peening  processes.  The  depth  of  compressive 
stress  is  directly  related  to  the  depth  of  plastically  affected  zone.  Ballard  et 
al.  [27,281  provided  an  empirical  expression  for  plastically  affected  depth  {Dp) 
and  surface  residual  stress  (Csurf)- 


f  VpVe^  \ 
\ye-Vp) 


^surf  — 


1+V 

1-V 


/  P-Ohel\ 

V  ^^HEL  ) 


4v^ 

71 


(4) 

(5) 


Here,  14,  Vp,  x,  P,  and  <5hel  indicate  the  elastic  wave  speed,  plastic  wave 
speed,  shear  stress,  shock  wave  pressure,  and  Hugoniot  Elastic  Limit  of  the 
material.  And  /d,  Zp,  a  indicate  pulse  duration,  plastic  surface  strain,  and  the 
side  of  square  impact.  Equations  4  and  5  indicate  plastically  affected  depth 
and  surface  residual  stress,  respectively.  The  LP  process  parameters  can  be 
optimized  to  generate  a  higher  volume  of  compressive  stress  with  constraints 
on  depth  of  compressive  stress  and  on  magnitude  of  tensile  stress. 


1.2.4  Initial  Developments 

Initially,  White  [291  identified  the  shock  wave  generation  ability  of  high-energy 
pulsed  lasers.  This  discovery  led  to  investigations  of  laser-generated  appli¬ 
cations  [301.  This  ability  was  explored  in  many  aspects  [31,321  including 
altering  material  properties  [331,  producing  vacancies  in  thin  vanadium  and 
nickel  foils  [341,  and  detonating  explosives  [351.  In  one  of  the  investigations 
by  Fairand  [361,  h  was  discovered  that  a  laser-generated  pressure  pulse  can 
alter  in-depth  material  properties  in  metals.  It  was  observed  that  the  yield 
strength  of  Aluminum  was  improved  by  laser-induced  shocks,  without  signifi¬ 
cantly  affecting  the  ultimate  tensile  strength  [361.  In  1974,  a  patent  was  issued 
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to  Malozzi  and  Fairland  [33]  for  LP’s  use  in  modification  of  material  prop¬ 
erties.  Upon  observing  that  a  higher  pressure  pulse  produces  a  higher  shock 
wave,  Anderhlom  [37]  established  that  higher  pressure  can  be  achieved  by 
applying  confined  ablation. 

After  these  inventions,  Clauer  et  al.  [1, 2, 38,  39]  and  Fairand  et  al.  [3-5,  36] 
performed  a  number  of  investigations  at  Battelle  Columbus  Laboratories,  Ohio 
(USA)  between  1968  and  1991,  exploring  the  potential  of  this  technology. 
Many  researchers,  such  as  Fox  [40],  Ford  et  al.  [41],  Ortiz  and  Penny  [42], 
and  Forget  [43],  also  contributed  to  the  development  of  the  process.  Fairand 
and  Clauer  [4,5]  delved  into  the  characteristic  of  the  laser-material  interaction 
and  the  material  response  to  the  induced  shock  waves.  In  this  research,  the 
peak  pressure  was  investigated  for  different  ablative  overlays  and  the  effects 
of  LP  on  welded  joints  was  reported.  After  these  investigations,  LP  was  ap¬ 
plied  to  improve  the  fatigue  life  of  fastener  holes  [6].  Fairand  and  Clauer  also 
performed  systematic  investigations  of  the  parameters  that  significantly  affect 
the  LP  process.  These  works  led  to  the  first  well-known  successful  indus¬ 
trial  use  of  LP  in  aircraft  gas  turbine  blades  [44]  for  increased  foreign-object 
damage  resistance.  Fox  [40]  investigated  the  effects  of  coating  on  the  target 
surface.  Ford  et  al.  [41]  showed  that  a  decrease  in  crack  growth  rate  can  be 
achieved  by  careful  selection  of  the  process  parameters  under  specific  condi¬ 
tions.  Ortiz  and  Penny  [42]  obtained  a  patent  for  developing  an  LP  system 
that  includes  a  foil  aligned  with  a  surface  of  the  workpiece  to  be  peened  and 
lasing  the  aligned  foil  surface.  This  system  improved  the  speed  of  peening  and 
accuracy  of  peening  control.  Forget  et  al.  [43]  studied  the  laser  beam  optical 
phenomena  and  deformation  mechanisms  and  proposed  an  analytical  model 
for  predicting  residual  stresses. 

1.2.5  Growth  Challenges 

The  success  of  initial  results  generated  interest  in  further  investigations  [38,39, 
41,45],  followed  by  industrial  applications  [42,46]  of  the  technology.  Many 
researches  showed  the  benefits  of  the  technology  and  produced  extensive  lit¬ 
erature.  After  the  initial  demonstration  of  the  technology  proved  it  to  be  ben¬ 
eficial  for  intended  material  properties,  the  number  of  LP  applications  grew 
rapidly  [47].  At  this  stage,  this  enabling  technology  faced  many  challenges. 
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The  first  challenge  was  to  develop  a  mechanism  that  could  deliver  a  powerful, 
controllable,  and  repeatable  pressure  pulse  at  a  high  frequency.  The  second 
challenge  was  to  develop  a  comprehensive  knowledge  of  the  process  and  a  nu¬ 
merical  tool  to  predict  the  generated  residual  stress  field.  The  understanding 
of  the  process  can  guide  the  research  in  the  field  and  could  assist  in  developing 
a  numerical  tool.  The  third  challenge  was  to  accurately  predict  and  document 
the  benefits  of  the  LP  process  on  material  strength.  Overcoming  these  chal¬ 
lenges  could  take  the  technology  to  the  next  level,  allow  an  optimal  use  of  the 
technology,  and  facilitate  reduction  in  the  cost  of  industrial  applications. 

This  section  presents  the  research  work  performed  to  overcome  the  first  chal¬ 
lenge.  The  pressure  generated  by  a  laser  pulse  provides  a  controllable  and 
reproducible  tool  for  practical  material  shock  processing  applications.  Re¬ 
search  in  the  confinement  regime  [37,40,43,48-511  and  laser  systems  [521 
provided  a  potential  solution  to  obtain  compatible  pressure  pulse  properties. 
Anderhlom  [371  observed  that  plasma  pressure  was  significantly  increased  by 
the  presence  of  confinement  (as  compared  to  no  confinement).  Fox  [401  fur¬ 
ther  developed  the  idea  of  confinement  regime  and  studied  the  effects  of  the 
coating  used  for  confinement.  The  investigations  of  confined  interaction  mode 
by  Clauer  et  al.  [481  and  Sano  et  al.  [501  indicated  improvement  in  the  pres¬ 
sure  pulse  magnitude  by  a  factor  of  5  (for  the  peak  pressure  level)  and  in  the 
pressure  pulse  duration  by  a  factor  of  3  compared  to  a  direct  ablation  mode. 
Hong  et  al.  [511  also  investigated  the  same  problem  by  researching  five  dif¬ 
ferent  kinds  of  confinement  overlays.  The  limitations  of  confinement  regime 
were  observed  in  the  form  of  dielectric  breakdown  phenomena  by  Foret  et 
al.  [43,491.  The  breakdown  seems  to  be  due  to  internal  breakdown  of  the 
transparent  material  such  that  the  energy  is  absorbed  primarily  inside  the  trans¬ 
parent  material  before  reaching  the  black  paint  -  target  interface.  Thus  dielec¬ 
tric  breakdown  of  the  confining  medium  limits  the  maximum  peak  pressure 
obtained  by  increasing  laser  energy  density.  Development  in  the  field  of  con¬ 
finement  regime,  along  with  inventions  in  laser  technology  [23,  241,  helped 
in  generating  an  increased  pressure  pulse.  The  laser  systems  that  meet  the 
process  requirements  are  not  readily  available  in  the  market,  but  they  can  be 
custom  built.  Further  developments  in  the  area  of  laser  systems  could  reduce 
the  setup  cost  for  the  LP  processing  facility. 

The  second  challenge  LP  faces  is  the  development  of  a  comprehensive  under- 
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standing  of  the  process.  The  knowledge  required  is  described  in  three  aspects: 
pressure  pulse  generation,  shock  wave  propagation  in  materials,  and  numeri¬ 
cal  simulation  of  the  process.  Within  pressure  generation,  Fabbro  et  al.  [53] 
described  the  confined  ablation  mode  with  a  three-step  process.  Using  this 
model  and  considering  the  plasma  to  be  a  perfect  gas,  the  scaling  law  of  the 
pressure  generation  (P)  can  be  estimated  by  the  following  relationship: 

where  Iq  is  the  incident  laser  power  density  in  (Giga  Watt / cm^),  P  is  the  pres¬ 
sure  in  (GPa),  Z  is  the  reduced  shock  impedance  between  the  target  and  the 
confining  medium  in  ohm,  and  a  is  the  efficiency  of  the  interaction.  The  values 
of  a  range  from  0.1  to  0.2.  More  recently,  Wu  and  Shin  [54,55]  developed  a 
mathematical  model  to  predict  pressure  pulse  properties  (magnitude,  duration, 
and  shape). 

For  shock  wave  propagation  in  materials.  Love  [56]  discussed  the  historical 
developments  in  the  field  of  stress  propagation  in  elastic  solids.  The  propaga¬ 
tion  theory  of  shock  and  plastic  waves  in  material  are  found  in  Kolsky  [25]. 
The  experimental  data  [9, 10]  available  in  the  literature  for  the  material  under 
investigation  is  up  to  10^/.?  strain-rate,  whereas  during  the  LP  process,  strain- 
rate  goes  up  to  lO^/j".  Further  advances  in  material  data  creation  will  assist 
modeling  and  simulation  activity. 

The  goal  of  performing  LP  is  to  improve  the  material  properties,  such  as  fa¬ 
tigue  and  fretting  strength.  Therefore,  a  formulation  is  needed  that  can  evaluate 
the  influence  of  LP  parameters  in  terms  of  resulting  material  properties.  This 
evaluation  can  be  performed  in  two  steps.  In  the  first  step,  the  residual  stress 
field  is  determined  from  the  LP  parameters.  In  the  second  step,  material  prop¬ 
erties  are  assessed  from  the  obtained  residual  stress  field.  This  evaluation  is 
the  third  challenge. 

In  the  literature,  the  analysis  and  simulations  are  limited  to  the  first  part  of 
the  process.  In  order  to  achieve  an  effective  residual  stress  field,  an  appropri¬ 
ate  selection  of  the  process  parameters  is  required.  To  select  a  set  of  optimal 
parameters,  an  accurate  residual  stress  prediction  methodology  is  essential. 
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Figure  6:  Laser  pulse  and  resulting  pressure  pulse  on  a  target  (Peyre  et.  al.  1996) 


Clauer  et  al.  [1,57],  Peyre  et  al.  [6,7],  and  Fabbro  et  al.  [53]  performed  numer¬ 
ous  experiments  to  determine  residual  stress  fields  for  different  LP  parameters 
and  materials,  and  showed  the  advantages  of  LP  in  improving  the  fatigue  life. 
In  his  research,  Peyre  et.  al.  obtained  a  pressure  pulse  shape  shown  in  Figure 
6.  Peyre  et  al.  also  performed  experiments  for  repeated  shots,  for  comparing 
residual  stress  field  induced  by  shot-peening  and  LP,  and  finally  for  compar¬ 
ing  fatigue  life  improvements  by  shot-peening  and  LP.  The  results  showed  that 
fatigue  life  improvement  through  laser  peening  is  higher  as  compare  to  shot 
peening.  Clauer  et  al.  [1]  obtained  the  pressure  pulse  properties  by  placing 
pressure  gauges  behind  the  surface  of  a  thin  specimen.  Peyre  et  al.  [58, 59] 
and  Ballard  [27]  obtained  the  pressure  pulse  properties  from  the  VISAR  (ve¬ 
locity  interferometer  system  for  any  refiector)  pressure  determination  system. 

Experiments  alone  are  not  sufficient  to  demonstrate,  expand,  and  advance  the 
needs  of  the  industry  because  LP  experiments  are  complicated,  potentially  te¬ 
dious,  and  expensive.  The  experimental  procedure  becomes  more  difficult  be¬ 
cause  of  the  use  of  destructive  stress  measurement  techniques  such  as  X-ray 
diffraction.  These  reasons  provide  a  motivation  to  develop  analytical  and  sim- 
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ulation  capabilities  of  the  process.  Clauer  et  al.  [57]  also  noted  that  having  a 
comprehensive  model  of  the  process  would  save  significant  time  and  cost  when 
developing  new  applications.  A  model  could  provide  guidance  and  what-if 
studies  in  modification  and  further  development  of  the  process.  In  addition, 
a  model  would  be  invaluable  in  enabling  close  process  control  in  production 
use  of  laser  peening.  The  ideal  scenario  would  be  to  have  a  system  that  needs 
a  geometric  model,  material  properties,  and  laser  system  configuration  and 
provides  the  engineer  with  optimal  process  parameters  to  achieve  the  desired 
residual  stress  field. 

In  modeling  the  LP  process,  the  prediction  of  the  resultant  residual  stress  field 
is  similar  to  the  analysis  of  shock  waves  in  materials.  An  analytical  model 
of  LP  was  first  developed  by  Ballard  et  al.  [28].  This  was  a  one-dimensional 
model  which  calculated  plastic  deformation  and  residual  stress  levels  as  a  func¬ 
tion  of  impact  pressure.  A  second  analytical  model  was  developed  by  Forget 
et.  al.  [43,49],  which  focused  on  the  surface  effects  induced  by  shock  waves. 
This  model  applies  to  circular  impacts  where  shear  waves  create  interference 
phenomena.  These  models  were  further  investigated  by  Peyre  et  al.  [6]  and 
Fabbro  et  al.  [53].  Peyre  at  al.  proposed  a  relation  between  the  peak  pressure 
pulse  magnitude  and  laser  power  density  (/q).  Fabbro  et  al.  [53]  discussed  the 
effects  of  laser  intensity,  target  materials  properties,  laser  pulse  duration,  and 
laser  wavelength  in  the  confinement  regime. 

These  analytical  efforts  are  based  on  elastic-perfectly  plastic  material  model 
assuming  uniaxial  strain  conditions.  This  assumption  is  fairly  accurate  for  the 
initial  shock  wave  response.  However,  under  many  conditions,  this  assumption 
may  not  hold  true  because  of  the  presence  of  non-uniaxial  strain  affecting  the 
resultant  residual  stress  field.  These  analytical  models  are  important  steps  in 
the  development  of  the  overall  LP  process,  however  not  sufficient  to  provide 
the  understanding  of  the  process  and  practical  applications. 

To  overcome  the  short  comings  of  analytical  modeling,  multiple  efforts  were 
involved  in  developing  FE  simulation  for  the  process.  After  extensive  experi¬ 
ments,  Nam  [9]  and  Noll  [10]  developed  their  own  FE  code  for  the  LP  simu¬ 
lations,  which  they  successfully  implemented  for  a  two-dimensional  process. 
Braisted  and  Brockman  [8]  performed  FE  analysis  of  one-sided  and  two-sided 
LP  for  an  axi-symmetric  2-D  case  using  ABAQUS  for  a  single  LP  treatment. 


21 


Ding  and  Ye  [11,60]  made  initial  attempts  in  the  3 -Dimensional  simulation; 
however,  they  performed  only  single  and  multiple  treatments  of  LP  at  the  same 
location  on  simple  geometries. 

In  the  process  of  developing  an  appropriate  simulation  methodology,  various 
researchers  have  used  different  parameters  depending  upon  the  experimental 
setup  at  hand.  The  key  differences  in  these  simulation  models  are  the  pressure 
pulse  shape  and  the  material  model.  Ding  and  Ye  [11,60]  used  a  perfecdy 
e]asto-p]astic  materia]  mode]  using  dynamic  yieid  strength  (HEL).  Zhang  et 
a].  [61]  investigated  the  simuiation  for  sma]]  impacts  with  the  Steinberg  mode], 
taking  pressure  effect  into  account  for  the  materia]  mode].  However  Peyre  et 
a].  [58,59]  found  constant  shock  yie]d  strength  HEL  at  different  LP  pressures. 
As  a  resuit,  they  used  the  Johnson-Cook  mode]  with  isotropic  hardening  which 
is  a  strain,  strain-rate,  and  temperature  dependent  mode].  In  this  research,  they 
also  presented  a  limited  investigation  of  the  influence  of  process  parameters 
such  as  pressure  pulse  amplitude  and  duration,  laser  spot  size,  and  sacriflcial 
overlay.  Recently,  Wu  and  Shin  [54,  55]  deveioped  a  mode],  this  inciudes 
anaiytica]  modehng  of  the  pressure  puise  from  a  iaser  beam,  and  iimited  PE 
simuiations  for  residua]  stress  prediction. 

The  numerica]  work  discussed  so  far  concentrated  on  deveioping  a  methodoi- 
ogy  for  PE  simuiation  of  the  process  for  a  iimited  number  of  the  LP  parameters 
and  simpie  geometries.  However,  to  determine  the  usefuiness  of  the  generated 
residua]  stress  fleid,  a  numerica]  process  is  needed  that  can  assess  materia] 
properties  with  or  without  the  presence  of  residua]  stresses.  The  deveiopment 
of  this  simuiation  research  can  reiate  the  process  parameters  to  the  fatigue  hfe. 
Therefore,  it  is  necessary  to  determine  the  improvements  in  the  fatigue  iife 
from  a  set  of  LP  parameters. 

1.2.6  LP  Advantages  and  Disadvantages 

This  section  iists  the  advantages  and  disadvantages  of  the  LP  process. 

•  It  is  wefl  known  that  the  cyciic  behavior  of  mechanicai  components  de¬ 
pends  on  the  metaiiurgicai,  mechanicai,  and  geometricai  surface  states 
of  the  materia]  invoived.  Many  experimentai  and  practicai  resuits  have 
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proved  that  fatigue  properties  are  improved  by  LP. 

•  A  critical  drawback  of  SP  is  that  residual  stresses  are  induced  by  the  con¬ 
tact  between  the  bombarding  ball  and  surface;  hence  residual  stresses  tend 
to  relax  quickly  under  repeated  loading.  Since  there  is  no  contact  during 
the  LP  treatment,  residual  stresses  tend  to  relax  more  slowly  than  SP. 

•  In  the  case  of  SP  or  LPB,  the  surface  of  the  ball  or  the  roller  tends  to 
degrade  with  prolonged  usage.  This  degradation  leads  to  a  rough  surface 
finish  on  the  target  component.  The  LP  process  does  not  use  any  physical 
tool  to  induce  residual  stress. 

•  As  mentioned  in  Section  1.2,  the  applicability  of  LP  to  complicated  ge¬ 
ometries  is  a  unique  advantage,  made  possible  because  a  laser  is  a  collec¬ 
tion  of  rays  that  can  reach  any  intricate  location  within  a  bulky  apparatus. 

•  The  cost  of  LP  is  relatively  high;  however,  the  cost  should  decrease  as  the 
advances  in  laser  technology  continue. 

•  For  a  successful  surface  enhancement  technique,  complete  control  of  the 
process  is  essential  to  generate  the  required  residual  stress  field.  Other¬ 
wise,  a  technique  may  carry  some  risk.  In  the  application  of  SP,  many 
probabilistic  parameters,  such  as  area  of  contact,  angle  of  contact,  and 
speed  of  strike,  are  involved.  These  make  the  process  difficult  to  control. 
LP  is  a  comparatively  controllable  and  repeatable  process. 

•  Although  LP  is  a  controllable  process,  it  tends  to  produce  non-uniform 
residual  stress  across  the  laser  spot.  This  non-uniformity  depends  upon 
the  metallurgical  properties  of  the  peened  component. 


1.3  Motivation 

In  the  past,  researchers  have  experimentally  investigated  the  effects  of  LP  pa¬ 
rameters  on  fatigue  [1^,6, 71  and  performed  a  limited  number  of  simulations 
on  simple  geometries  [8-10,54,55,59,601.  However,  experimentally  monitor¬ 
ing  the  dynamic,  intricate  relationship  of  peened  material  is  challenging.  With 
an  increasing  number  of  applications  [62,  631  to  complex  geometries  in  the 
nuclear  industry,  to  aircraft  turbine  blades,  and  to  medical  implants,  these  lim¬ 
ited  experimental  and  simulations  capabilities  are  not  sufficient  for  effective 
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LP  process  design.  A  comprehensive  procedure  is  required  that  can  perform 
simulations  of  multiple  treatments  of  LP  at  the  same  location  and  sequential 
LP  at  multiple  locations,  different  overlapping  configurations  of  LP  locations, 
and  application  on  complex  geometries.  With  increased  computer  speed  as 
well  as  increased  sophistication  in  non-linear  finite  element  analysis  (FEA) 
softwares,  it  is  now  possible  to  develop  a  model-based  design  optimization 
approach  for  an  effective  LP  application. 

For  a  typical  job,  the  component  is  peened  using  a  few  (usually  4)  parameter 
settings.  These  parameter  settings  for  the  experiments  are  selected  based  on 
experience.  After  peening,  all  the  components  are  subject  to  service  loads  to 
which  components  are  subject  to  during  their  intended  application.  Based  on 
the  improvement  achieved  in  different  components,  the  actual  process  parame¬ 
ters  are  selected  for  the  job.  Depending  upon  the  resources  available,  the  above 
process  can  be  repeated  after  refining  the  parameters. 

For  the  LP  industry  to  move  forward  with  applications,  it  is  necessary  to  take 
the  simulation  to  the  next  level  from  simple  simulation  case  studies.  There  are 
several  parameters  that  demand  exploration  of  the  design  space  to  optimize 
the  residual  stress  field.  For  example,  the  spot  size,  spot  shape,  pulse  dura¬ 
tion,  pulse  shape,  pulse  intensity,  number  of  treatments,  sequence  of  locations, 
amount  of  overlap,  and  overall  LP  location  layout  can  each  make  a  significant 
difference  in  component  life.  Some  of  these  variables  are  discrete  and  oth¬ 
ers  are  continuous,  creating  a  mixed  variable  problem  for  optimization.  Along 
with  achieving  the  desired  residual  stress  field,  reduction  in  unfavorable  tensile 
region  (beneath  compressive  region)  is  also  a  goal.  Although  it  is  not  possible 
to  eliminate  the  entire  tensile  region,  it  can  be  constrained  to  be  less  trou¬ 
bling.  Solving  a  practical  mixed  design  variable  optimization  problem  will  be 
a  challenge.  The  problem  formulation  is  complicated  by  the  time-dependent 
nature  of  plastic  deformation  during  the  LP  process.  In  addition,  the  compu¬ 
tational  effort  needed  with  LP  simulations  is  significant  because  of  the  elastic 
and  plastic  analysis  involved  with  small  time  increments.  Therefore,  an  ef¬ 
ficient  simulation-based  methodology  is  needed  to  effectively  design  the  LP 
process.  Once  the  methodology  is  fully  developed,  it  can  be  used  to  optimize 
complicated  applications  for  direct  cost  savings. 
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1.4  Report  Organization 


This  report  is  organized  in  the  following  order: 

•  Section  1  describes  the  step-by-step  LP  process  in  detail.  This  includes 
brief  description  of  laser  generation,  overlay  application,  shock  wave  prop¬ 
agation,  and  residual  stress  generation.  The  physics  behind  the  generation 
of  residual  stress  is  discussed  followed  by  a  brief  outline  of  the  processes 
development  and  growth. 

•  Section  2  shows  the  different  modeling  parameters  of  the  LP  process,  and 
available  simulation  methodology  from  the  literature.  A  modified  method¬ 
ology  is  presented,  which  is  then  used  model  a  sequential,  multiple-location 
LP  treatment,  a  parametric  plate  model,  and  a  coupon  model.  Once  de¬ 
veloped,  this  methodology  can  be  used  to  perform  LP  simulations  of  a 
representative  realistic  component. 

•  Section  3  presents  the  background  of  optimization  and  the  issues  faced 
when  solving  the  mixed- variable  design  optimization  problem  of  the  LP 
process.  Gradient-  and  non-gradient-based  optimization  methods  are  in¬ 
troduced.  The  simulation  developed  in  section  3  is  employed  to  perform 
parametric  studies  of  the  LP  parameters. 

•  Section  4  discuss  the  idea  that  one  location  laser  peening  is  a  multimodal 
mixed- variable  optimization  problem.  To  effectively  solve  the  problem  a 
multimodal  mixed- variable  optimization  method  is  proposed.  This  method 
is  tested  on  multimodal  problems  and  successfully  implemented  to  solve 
the  LP  problem. 

•  Section  5  presents  a  strategy  to  optimize  multi-location  laser  peening.  The 
proposed  strategy  employs  multi-fidelity  and  surrogate  models  to  address 
the  issue  of  computational  cost. 

•  Section  6  concludes  the  report  with  the  list  of  contributions,  summary,  and 
future  work. 
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1.5  Chapter  Summary 


This  section  presents  surface  enhancement  techniques,  justification  of  LP,  and 
a  brief  overview  of  LP.  LP  process  details  are  discussed,  such  as  the  opaque 
and  transparent  overlays,  laser  generation,  shock  wave  propagation,  and  resid¬ 
ual  stress  generation.  A  brief  history  of  the  processes  development  and  growth 
is  reviewed.  In  the  overall  research,  this  section  serves  to  demonstrate  a  com¬ 
prehensive  understanding  of  the  process,  the  physics  involved,  and  the  moti¬ 
vation  behind  this  research. 


26 


2.0  FINITE  ELEMENT  SIMULATION 


Utilization  of  the  promising  surface  enhancement  technology  will  be  difficult 
with  only  the  knowledge  gained  through  experiments.  An  increasing  number 
of  applications  as  well  as  a  quest  to  understand  the  process,  demand  a  com¬ 
putational  exploration  of  the  process.  The  physical  and  historical  informa¬ 
tion  from  sections  1  and  2  will  help  in  developing  a  well-organized  simulation 
methodology. 

Before  beginning  the  simulation  process,  it  is  vital  to  research  the  different 
parameters  required  for  FE  simulation.  Apart  from  the  simulation  procedure, 
there  are  four  major  modeling  parameters.  These  are  the  laser  spot  size  and 
shape,  pressure  pulse  shape  and  duration,  geometric  model  and  meshing,  and 
material  model.  These  parameters  are  determined  based  on  published  experi¬ 
mental  and  simulation  results.  This  section  is  divided  into  three  major  sections. 
In  the  first  section,  a  brief  description  of  the  simulation  parameters  and  criteria 
for  determining  them  are  presented.  The  second  section  provides  an  expla¬ 
nation  of  the  simulation  procedure  used.  In  the  third  section,  the  simulation 
results  of  single  treatment  and  multiple  treatments  of  LP  at  the  same  location 
are  compared  with  published  experimental  results.  Once  the  simulation  pro¬ 
cess  is  validated  with  published  experimental  results,  geometric  shapes  such 
as  a  rectangular  box  are  analyzed  for  sequential  LP  at  multiple  locations.  Sim¬ 
ulation  results  for  two  configurations  (0%  and  50%  overlap)  of  LP  at  multiple 
locations  are  compared  with  each  other. 

This  is  followed  by  a  development  of  a  cylindrical  model,  parametric  plate, 
and  rectangular  bending  coupon  models.  The  cylindrical  model  is  employed 
to  study  the  effect  of  peening  sequence.  This  is  also  used  to  investigate  effects 
of  a  peening  shot  on  near-by,  previously  peened  area.  The  parametric  plate 
is  employed  to  optimally  design  a  residual  stress  field  for  a  fiat  surface.  The 
rectangular  coupon  model  can  be  used  not  only  for  peening  simulations,  but 
also  to  perform  fatigue  life  calculations. 
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2.1  Simulation  Parameters 


There  are  four  major  modeling  parameters  in  LP  FE  analysis.  These  are  the 
laser  spot  shape  and  size,  the  pressure  pulse  shape  and  duration,  geometric 
modeling  and  meshing,  and  the  material  model.  Apart  from  the  mesh  conver¬ 
gence  study,  these  parameters  are  selected-based  on  the  literature.  Details  of 
these  parameters  are  discussed  in  the  following  sections. 

2.1.1  Laser  Spot  Shape  &  Size 

Laser  spot  shape  and  size  are  vital  simulation  parameters.  These  shapes  affect 
the  magnitude  and  depth  of  residual  stress  and  can  determine  the  total  cost. 
Two  commonly  used  spot  shapes  are  circular  and  rectangular.  Other  shapes, 
such  as  elliptical  (and  its  variations)  and  rectangular  can  be  generated  from 
circular  and  square  laser  system  by  changing  the  striking  angle  of  laser  beam. 
The  accurate  effects  of  spot  shape  changes  are  not  known  and  further  explo¬ 
rations  of  the  parameter  are  needed.  As  far  as  the  size  of  the  spot  is  concerned, 
the  shock  wave  attenuation  rate  for  small  diameter  (0.5-1  mm)  spots  is  higher 
compared  to  large  diameter  (4-5  mm)  spots.  This  is  due  to  two-dimensional  at¬ 
tenuation  of  the  shock  wave  in  the  case  of  small  impacts  that  reduce  the  depth 
of  the  plastically  affected  region  and  of  the  compressive  residual  stress  field. 

The  selection  of  the  parameter  is  dependent  on  the  laser  generation  mech¬ 
anism,  material  behavior,  component  geometry,  and  job  requirements.  Al¬ 
though  it  is  possible  to  change  the  spot  shape  and  size  during  the  process,  one 
size  and  shape  is  used  for  the  entire  component,  due  to  the  effort  and  time 
involved  in  changing  these  parameters.  At  the  same  time  further  numerical  ex¬ 
ploration  of  various  shapes  is  necessary  to  estimate  their  benefits  before  they 
can  be  used  in  practical  applications. 

2.1.2  Pressure  Pulse  Magnitude  and  Shape 

Residual  stress  magnitude  and  depth  are  the  most  sensitive  to  pressure  pulse 
magnitude  and  shape  parameters.  The  magnitude  of  the  pressure  pulse  can  be 
controlled  by  the  intensity  of  the  laser  beam.  The  relation  between  laser  inten- 
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sity  and  the  peak  pressure  magnitude  in  confined  ablation  mode,  as  estimated 
by  Fabbro  et  al.  [53]  and  Peyre  et  al.  [6,58]  with  certain  assumptions,  is  given 
in  Equation  6.  The  pressure  pulse  intensity,  shape,  and  duration  depend  upon 
the  properties  of  the  laser  beam.  In  the  past,  experiments  have  been  performed 
to  determine  laser  beam  and  pressure  pulse  properties.  Figure  6  shows  typical 
shapes  of  the  laser  beam  pulse  and  pressure  pulse.  Depending  upon  various 
parameters,  the  laser  beam  midspan  varies  from  5-30  nanoseconds  (ns),  and 
the  pressure  pulse  midspan  is  approximately  three  times  that  of  the  laser  beam 
midspan. 

The  temporal  distribution  of  the  pressure  pulse  on  the  surface  of  the  target  ma¬ 
terial  can  be  experimentally  measured  by  multiple  methods.  Clauer  et  al.  [1] 
obtained  the  pressure  pulse  properties  by  placing  pressure  gauges  behind  the 
surface  of  a  thin  specimen.  Peyre  et  al.  [58,  59]  and  Ballard  [27]  obtained 
the  pressure  pulse  properties  from  the  VISAR  pressure  determination  system. 
VISAR  developed  by  Barker  and  Hollanbach  [64]  is  an  optical-based  system 
that  utilizes  Doppler  interferometry  technique  to  measure  the  time-history  of 
the  motion  of  a  surface.  The  laser  light  is  focused  to  a  point  on  a  target  of 
interest  and  the  reflected  light  is  collected,  routed  through  an  unequal  leg  in¬ 
terferometer,  and  converted  to  electronic  information.  This  information  is  then 
analyzed  for  the  amount  of  Doppler  shift  during  a  given  time  [65].  VISAR  can 
be  applied  to  a  wide  range  of  experimental  conditions  and  provides  accurate 
measurement  with  1  to  2  ns  time  resolution. 

Recently,  Wu  and  Shin  [54, 55]  developed  an  analytical  relation  between  the 
laser  parameters  and  the  LP  pressure  pulse.  Ding  and  Ye  [60]  assume  that 
shape  of  the  pressure  pulse  is  like  a  Gaussian  distribution  and  approximate  it 
using  a  triangular  pressure  pulse.  During  various  experiments,  it  was  observed 
that  a  steep  rise  in  pressure  tends  to  generate  a  better  residual  stress  field  in 
terms  of  magnitude  and  depth  as  compared  to  a  Gaussian  pulse.  Nam  [9]  and 
Noll  [10]  use  a  shape  as  shown  in  Figure  7.  The  figure  shows  the  ratio  of 
pressure  at  time  t  to  peak  pressure  of  the  pulse.  This  pressure  pulse  is  used  for 
experimental  validation  purpose.  The  pressure  pulse  used  in  this  research  for 
various  parametric  studies  and  optimization  is  little  different  than  this  shape. 
Overall  pressure  pulse  features  are  kept  similar.  But  minor  modifications  are 
made  after  investigating  the  pressure  pulse  shapes  available  in  the  literature. 
The  difference  between  the  selected  shape  and  this  shape  is  the  the  selected 
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Figure  7;  Temporal  loading  profile  of  pressure  pulse  (Nam,  2002) 
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Laset  Spot 


Figure  8:  Schematic  of  2-D  axi-symmetric  FE  model 


shape  is  little  smoother  profile  than  this  pulse. 

2.1.3  Geometric  Modeling  and  Meshing 

Normally,  the  LP-affected  zone  is  small  in  size  compared  to  the  component 
size.  At  the  same  time,  LP  is  a  high  speed  dynamic  process,  so  a  fine  mesh 
is  required  to  model  stress  waves  and  the  interactions  between  them.  A  small 
area  of  the  component  that  is  affected  by  LP,  called  the  control  volume,  is 
modeled  using  finer  elements,  and  the  area  surrounding  this  control  volume 
is  modeled  using  infinite  elements.  Infinite  elements  are  normally  used  to 
define  unbounded  domain  and  to  provide  a  model  with  a  ‘quiet’  boundary. 
Infinite  elements  are  used  in  a  model  in  which  the  region  of  interest  is  small 
in  size  compared  to  its  surroundings,  such  as  in  LP  modeling.  Figures  8  and 
9  show  a  two-  and  three-dimensional  representative  models  of  a  component. 
In  both  cases,  the  laser  spot  is  a  circle  of  5  mm  in  diameter.  The  experimental 
results  [91  show  spatial  variation  of  the  pressure  across  a  peening  spot.  To 
incorporate  spatial  variation,  as  shown  in  Figure  10,  the  spot  is  divided  into 
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Table  1;  Basic  material  data  for  Ti-6A1-4V  (Nam,  2002) 


Material  properties,  units 

Value 

Young’s  modulus,  E  [GPal 

113.8 

Poisson’s  ratio,  V 

0.342 

Density,  p  [kg/m^] 

4500 

Initial  yield  stress,  Tq  [MPal 

924 

Hugoniot  Elastic  Limit  [MPal 

2800 

five  sections.  The  first  section  is  a  small  circle,  (blue  in  Figure  10)  from  which 
the  remaining  sections  or  bands  (yellow,  light  green,  red,  and  brown),  radiate 
outward  in  concentric  circles.  The  center  most  region  has  minimum  pressure 
(0.8  time  the  peak  pressure)  and  the  outer  most  has  maximum  pressure  (1.0 
time  of  the  peak  pressure). 

2.1.4  Material  Modeling 

LP  generates  strain-rates  exceeding  10^  s~^  within  the  target  material.  The 
material  model  plays  a  crucial  role  in  accurately  simulating  a  process  with  such 
a  high  strain-rate.  In  the  literature,  researchers  [11,58,59,611  have  used  many 
material  models.  In  this  research,  two  material  models  were  explored  for  use 
in  LP  simulation.  The  first  model  uses  perfectly  elastic-plastic  material  [111 
properties.  In  this  model.  Young’s  modulus  and  dynamic  yield  strength  are 
used  to  define  material  properties.  The  dynamic  yield  strength  depends  upon 
the  Hugoniot  Elastic  Limit  of  the  material,  and  is  calculated  using  Equation 
1.  The  material  properties  of  Ti-6A1-4V  used  for  all  simulations  are  shown  in 
Table  1  [101.  In  the  second  model,  to  model  the  stress-strain  dependence  at 
high  strain  rates,  the  Johnson-Cook  model  [59,66-681  is  used.  In  the  Johnson- 
Cook  model,  Ogq  is  expressed  as 


G,^  =  (A  +  5£”^) 


1+C.Ln 


(7) 


where  Oeq  is  the  effective  stress,  A,B,C  and  m  are  material  constants,  n  is 
the  work  hardening  exponent,  is  the  normalized  effective  plastic  strain  rate. 
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Figure  9:  Schematic  FEA  model  of  a  quarter  section 


Figure  10;  Schematic  of  discretized  spot  for  spatial  pressure  variation 
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Table  2:  Parameters  for  Johnson-Cook  material  model  of  Ti-6A1-4V  (Lesuer,  2000) 


A  (MPa) 

B  (MPa) 

n 

C 

m 

Ti-6A1-4V 

1098 

1092 

0.93 

0.014 

1.1 

and  the  quantity  T*  is  defined  as 


T 


* 


T-Tq  ' 

Tmelt  -  To, 


(8) 


where  Tf^eit  is  the  melting  temperature.  The  strength  of  the  material  is  a  func¬ 
tion  of  strain,  strain-rate,  and  temperature.  The  values  of  A,  B,C,n,  and  m  are 
determined  from  an  empirical  fit  of  experimental  flow-stress  data.  The  opaque 
overlay  prevents  the  material  from  generating  high  temperatures  during  the 
peening  process.  Due  to  this  coating,  the  effect  of  temperature  during  peen- 
ing  is  limited.  Therefore,  the  effects  of  temperature  are  not  considered  in  the 
material  model. 


2.2  Simulation  Procedure 

Because  of  the  complexity  of  elastic  and  plastic  deformation  and  shock  wave 
propagation,  it  is  necessary  to  employ  an  accurate  procedure  with  manageable 
computational  requirements.  As  discussed  in  section  1,  very  high  magnitude 
pressure  pulse  of  the  LP  process  lasts  for  very  small  period  of  time.  However, 
it  takes  a  much  longer  time  for  the  material  reach  an  equilibrium.  The  time 
scale  of  pressure  duration  is  in  nano  second  (ns).  The  time  scale  to  reach  equi¬ 
librium  is  in  milliseconds  (ms).  In  cases  of  shock  loading  such  as  in  LP,  equi¬ 
librium  can  not  be  achieved  quickly  because  of  the  propagation  and  interaction 
of  stress  waves.  The  residual  stresses  are  caused  by  the  plastic  deformation, 
and  the  plastic  process  last  for  a  very  short  time.  It  is  essential  to  obtain  the 
equilibrium  state  for  accurate  results. 

An  appropriate  simulation  procedure  is  needed.  The  explicit  algorithm  is  suit¬ 
able  for  high  speed  process  [691  such  as  initial  phase  of  the  LP  process,  but 
if  this  algorithm  is  used  for  the  entire  process  till  equilibrium,  it  is  not  com¬ 
putationally  efficient.  The  reason  for  the  inefficiency  is  that  this  algorithm  is 
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Figure  1 1 ;  Flow-chart  of  simulation  process 


not  designed  for  use  in  slow  processes  such  as  later  part  of  the  LP  process 
and  it  also  integrates  over  time  in  small  time  steps  for  algorithm  stability.  The 
implicit  algorithm  [691  of  FE  analysis  is  suitable  for  the  later  part,  but  not 
the  initial  part,  of  the  process.  Therefore,  a  simulation  procedure  is  used  that 
takes  advantage  of  both  processes  and  obtains  the  equilibrium  state  within  a 
manageable  computational  time. 

FE  simulation  of  the  LP  process  is  divided  into  two  phases  [81.  The  flow  chart 
of  the  process  is  shown  in  Figure  11.  The  flrst  is  modeling  the  initial  high  speed 
transient  high  pressure  loading  in  ABAQUS/  Explicit,  which  is  speciflcally 
designed  for  such  processes.  The  second  phase,  which  is  a  comparatively 
slower  physical  process,  is  to  determine  the  flnal  equilibrium  residual  stress 
held.  This  phase  of  the  analysis  is  performed  in  ABAQUS/Standard  (implicit 
algorithm). 
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2.2.1  Explicit  and  Implicit  Algorithms 


To  analyze  the  dynamic  response  of  a  material,  caused  by  high  speed  and 
high  intensity  laser  shock,  the  explicit  FE  method  is  used.  The  governing 
equations  for  general  explicit  algorithms  can  be  derived  by  requiring  the  work 
of  external  forces  to  be  absorbed  by  the  work  of  internal  and  inertia  forces  for 
any  small,  kinematically  admissible  motion,  which  satisfies  both  compatibility 
and  essential  boundary  conditions  [70,711: 


pM/6M/<if2+  /  X)puihuid£L 


Q 


Q 


foij5DijdO.  (9) 

Q. 


where  p  is  the  mass  density,  il  and  li  are  the  nodal  acceleration  and  velocity, 
respectively,  t)  is  the  damping  coefficient,  5u  is  the  virtual  displacement,  //  is 
the  body  force  density,  Tr  is  the  boundary  force  applied  on  the  boundary  F,  Oij 
is  the  Cauchy  stress  tensors,  and  Dij  =  is  the  deformation-rate 

tensor.  In  Eq.  9,  the  five  terms  represent  inertia,  damping,  body  forces,  bound¬ 
ary  forces,  and  strain  energy  of  the  system  under  consideration,  respectively. 
The  standard  FE  procedure  can  be  used  to  obtain  a  discretized  equation: 


[M]{U}  +  [C]{U}  +  [K]{U}  =  (10) 

where  [M]  is  the  diagonal  mass  matrix,  {U}  is  the  nodal  acceleration  vector, 
{U}  is  the  velocity  vector,  [C]  is  the  diagonal  damping  matrix,  is  the 

internal  element  force  vector,  and  is  the  external  and  body  force  vectors. 

In  the  integration  methods,  response  history  is  calculated  using  step-by-step 
integration  in  time.  The  response  is  evaluated  at  time  instants  At,  2 At,  ...,  nAt 
and  so  on.  At  the  time  step,  the  equation  of  motion  is  [711 


[M]  {U}„  +  [C]  {U}„  +  [K]  {U}„  =  {F^^%  (11) 


Discretization  in  time  is  achieved  by  finite  difference  approximation  of  time 
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derivatives.  There  are  many  ways  to  perform  the  discretization  and  solve  the 
equation.  In  this  discussion,  [M]  is  assumed  positive  definite  and  [K]  is  positive 
semidefinite.  Structure  is  allowed  to  have  rigid-body  motion  as  part  of  its 
response. 

In  the  time  integration  conditions  at  time  step  (t  +  1 )  are  calculated  from  the 
equation  of  motion,  a  difference  expression,  and  known  conditions  at  one  or 
more  preceding  time  steps.  Algorithms  can  be  classified  as  explicit  and  implicit 
depending  upon  the  information  used  in  calculating  conditions  at  (t  +  1 ) .  An 
explicit  algorithm  uses  an  expression  of  the  general  form  [711 


{U}„+1  =  /({U}„,  {U}„,  {U}„,  {U}„_1, . . .)  (12) 

that  contains  only  historical  information  on  its  right  hand  side.  This  equation 
is  combined  with  the  governing  Eq.  11  at  step  n  to  solve  the  problem.  An 
implicit  algorithm  uses  the  following  general  form 


{U}„+1  =  /({U}„+1,  {U}„,  {U}„,  {U}„,  . . .)  (13) 

that  is  combined  with  the  governing  Eq.  1 1  at  step  («+ 1)  to  solve  the  problem. 
In  the  explicit  algorithm,  only  the  historical  information  is  used  to  calculate 
dynamic  quantities  at  (t  +  At).  The  central  difference  method  is  used  in  the 
explicit  algorithm.  By  manipulating  the  Taylor  series  expansion  of  {U}„+i 
and  about  time 

{U}»  =  -  {U}„_i)  (14) 


It  can  be  written  as: 


{U}„+i  =  {U}„_i  +  2At{U}„ 


The  {U}„+i  can  be  used  to  calculate  {U}«: 


(15) 
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(16) 


{U„}  =  ^({U}„+,  -2{U}„  +  {U}„_,) 

The  integration  rule  in  the  explicit  algorithm  is  simple,  this  is  a  reason  behind 
the  computational  efficiency  of  the  explicit  dynamics  procedure.  The  equation 
derived  from  the  central  difference  scheme  used  in  the  algorithm  is  condition¬ 
ally  stable  [711.  As  a  result  stable  time  step  is  an  important  factor  and  the 
explicit  FEA  integrates  through  time  by  using  many  small  time  increments. 
The  stability  limit  for  the  time  step  operator  depends  upon  the  highest  natural 
frequency.  Since  it  is  hard  to  find  the  highest  natural  frequency,  a  conserva¬ 
tive  estimate  of  the  stable  time  increment  is  given  in  Equation  17  [69,711  in 
terms  of  minimum  element  size  (Lg)  and  the  current  effective  dilatational  wave 
speed  of  the  material  (Cj).  The  point  to  note  is  that  if  the  minimum  element 
size  is  reduced,  then  Atstabie  decreases  while  the  time  required  to  perform  the 
integration,  and  ultimately  analysis  time,  increases. 


At  St  able  (17) 

Due  to  this  limit  (17),  a  very  small  time  step  for  integration  is  used  for  ac¬ 
curate  results,  which  leads  to  large  computational  cost.  However,  the  implicit 
algorithm  removes  this  upper  bound  on  the  time  step  for  integration  by  solving 
for  the  dynamic  quantities  at  [t  +  At)  based  not  only  on  values  at  t,  but  also 
on  these  same  quantities  at  {t  ^  At).  Based  on  the  special  case  of  Newmark 
method  [721  the  relations  for  implicit  algorithm  as  give  by  [711  are 

{U}„+i  =  {U}„  +  A((7{U}„+i  +  (1-7){U}„)  (18) 


{U}„+i  =  {U}„  +  A<{U}„  +  [2p{U}„+i  +  (1  -  2|5){U}„]  (19) 


where  y  and  P  control  characteristics  of  the  algorithm  such  as  accuracy,  numer¬ 
ical  stability,  and  the  amount  of  algorithm  damping.  Respectively,  the  average 
acceleration  and  linear  acceleration  methods  [711  are  given  by  y^^^  P  =  ^  and 
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by  P  =  g.  By  solving  Eqs.  18  and  19  for  {U„+i},  then  further  substitu¬ 
tion  [711  can  provide  the  implicit  solution. 

It  is  difficult  to  solve  implicit  equation  in  the  algorithms,  but  the  time  step  can 
be  two  orders  of  magnitude  larger  than  the  explicit  algorithm.  Thus,  the  rela¬ 
tive  computational  cost  of  the  two  algorithms  depends  on  the  stability  limit  of 
the  explicit  scheme,  the  ease  with  which  the  nonlinear  equations  can  be  solved 
for  the  implicit  operator,  the  relative  size  of  time  increments  that  can  provide 
acceptable  accuracy  with  the  implicit  scheme  compared  to  the  stability  limit 
of  the  explicit  scheme,  and  the  size  of  the  model.  The  equilibrium  equation 
for  the  implicit  algorithm  can  be  obtained  from  Refs.  [69,711.  A  brief  de¬ 
scription  is  provided  in  this  section.  The  description  is  for  academic  purpose 
only.  These  algorithms  are  in-built  in  the  FEA  software.  The  user  only  need 
to  selected  the  appropriate  option  to  use  implicit  or  explicit  algorithm. 

Because  of  the  two  steps  involved  in  the  simulation  procedure,  it  is  necessary 
for  accuracy  and  efficiency  to  determine  the  analysis  duration  of  each  step.  The 
simulation  duration  would  be  much  longer  than  the  pressure  pulse  duration. 
Ding  and  Ye  [601  suggested  that  the  solution  time  should  be  two  orders  of 
magnitude  higher  (  3000)  than  the  pressure  pulse  duration  (  30).  In  this  work, 
a  large  number  of  analyses  is  performed  to  select  the  duration,  and  an  analysis 
time  of  2500  to  3000  ns  for  the  first  phase  and  10  —  50  millisecond  (ms)  for 
the  second  phase  is  used. 

2.2.2  Ohio  Supercomputing  Center  for  Simulations 

Due  to  the  computational  cost  involved  in  a  LP  simulation,  most  of  the  finite 
element  simulations  are  performed  on  the  Glenn  cluster  at  the  Ohio  Super¬ 
computing  Center  (OSC).  The  OSC  is  located  in  Columbus,  Ohio.  The  details 
of  the  Glenn  and  other  clusters  at  OSC  can  be  found  at  the  institute  website 
(www.osc.edu).  In  brief  [731,  there  are  more  than  7000  nodes,  each  node  with 
dual  or  quad  core,  2.6  GHz  speed,  and  8GB  or  higher  (more  than  6000  nodes 
with  8GB,  70  nodes  with  16GB,  16  nodes  with  32GB,  2  nodes  with  64GB, 
76  nodes  with  218GB,  and  10  nodes  with  1.8TB  RAM)  RAM.  The  total  CPU 
computational  time  at  OSC  involves  two  parts.  The  first  part  is  the  waiting 
period  and  the  second  part  is  the  computation  part.  The  waiting  time  is  a  of 
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function  many  parameters  such  as  load  on  server,  requested  CPU  time,  re¬ 
quested  RAM,  and  cluster  maintenance  activities.  The  CPU  time  mentioned 
in  the  report  does  not  include  the  waiting  time. 

2.3  Two-Dimensional  Simulation 

The  2D  simulation  is  an  efficient  way  to  simulate  one  shot  with  the  circular 
spot  at  the  same  location  because  of  less  computational  requirements.  Figure 
8  shows  an  axisyymmetric  FE  model  of  a  component.  The  model  consists  of  a 
5x5  (mm^)  rectangular  area  of  finer  elements  surrounded  by  5  mm  of  infinite 
elements  on  each  of  the  two  sides.  The  FE  model  represents  a  cylindrical  part 
of  10  mm  diameter  with  a  laser  spot  of  5  mm  diameter.  Based  on  the  mesh  con¬ 
vergence  a  200  X  200  mesh  is  used  for  the  5  mm  square.  This  model  contains 
40000  finer  elements  (four-noded  CAX4R  elements  in  ABAQUS)  and  400  in¬ 
finite  elements  (four-noded  CINAX4  elements).  This  model  is  developed  for 
use  in  solving  the  optimization  problem  of  designing  a  required  residual  stress 
field. 

2.4  Three-Dimensional  Simulation 

Using  2D  FE  simulations  is  a  computationally  efficient  way  to  predict  the 
residual  stress  field  for  one  shot  of  LP  at  the  same  location.  However,  vari¬ 
ations  such  as  sequential  shots  at  multiple  locations  or  rectangular  or  ellip¬ 
tical  spot  shapes  are  not  possible  to  simulate  in  2D.  These  variations  require 
three-dimensional  FE  simulation  of  the  process.  Figure  9  shows  a  quarter  FEA 
model  of  a  component.  The  quarter  model  consists  of  a  5-mm  cube  of  finer  el¬ 
ements  surrounded  by  5  mm  of  infinite  elements  on  each  of  the  three  sides. 
The  laser  spot  is  a  circle,  5  mm  in  diameter. 

2.4.1  Mesh  Convergence 

To  determine  a  suitable  mesh  size  for  the  simulation,  a  mesh  size  convergence 
study  is  performed  on  a  quarter  model  (Figure  9).  The  mesh  convergence  study 
is  divided  into  two  parts.  In  the  first  part,  four  uniform  meshes  are  investigated. 
In  the  second  part,  based  on  the  results  of  the  first  part,  three  biased  meshes 
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are  investigated.  The  first  part  is  to  determine  the  converged  mesh,  and  the 
second  part  is  to  bring  the  computational  cost  within  acceptable  limits  for  the 
research  purpose.  In  all  the  models,  control  volume  is  modeled  using  eight- 
noded  C3D8R  elements  and  the  surrounding  region  is  modeled  using  CIN3D8 
infinite  elements.  In  mesh  convergence,  40,  80,  and  160  elements  are  used  to 
mesh  a  5-mm-long  side  of  the  control  volume.  The  element  size  for  each  mesh 
is  0.125  mm,  0.0625  mm,  and  0.03125  mm,  respectively. 

The  mesh  size  with  320  elements  (0.015625  mm)  was  also  attempted.  It  is 
infeasible  to  use  this  mesh  due  to  the  computational  space  and  memory  re¬ 
quirements.  Deu  to  this,  mesh  size  with  200  element  was  attempted.  Even 
this  mesh  will  take  approximately  10  days  of  CPU  time.  This  time  does  not 
include  waiting  time  of  around  3-20  days.  The  waiting  time  is  high  because  of 
a  limited  number  of  high  memory  nodes  at  OSC. 

A  comparison  of  the  residual  stress  profile  from  various  mesh  sizes  is  shown  in 
Figure  12.  The  results  show  residual  stress  magnitude  along  T-axis  and  depth 
along  X-axis.  The  figure  compares  the  results  from  uniform  and  nonuniform 
mesh  models.  The  residual  stress  profile  can  be  divided  into  two  regions;  0.0 
to  0.2  mm  and  0.2  mm  to  5  mm.  The  residual  stress  profiles  between  the  depth 
of  0.2  to  5  mm  are  matching  with  each  other  for  most  of  the  mesh  sizes.  In  fact, 
all  uniform  mesh  model  results  are  on  top  of  each  other.  The  residual  stress 
profiles  between  the  depth  of  0.0  to  0.2  mm  are  different  from  each  other  for 
the  considered  mesh  sizes.  The  maximum  difference  between  mesh  sizes  80 
(80  X  80  X  80)  and  160  is  16%.  This  is  the  maximum  difference  and  is  large  but 
it  is  only  in  very  small  region.  The  results  are  matching  with  each  other  in  most 
of  the  depth.  Therefore,  the  160  mesh  can  be  assumed  to  be  a  converged  mesh, 
CPU  time  required  for  this  simulation  is  prohibitively  expensive  (7.5  days). 
The  use  of  this  mesh  will  make  it  difficult  to  perform  full  3Z)  simulations  for 
multiple  peening  locations. 

Based  on  the  time  taken  by  the  finest  mesh,  biased  mesh  convergence  study  is 
also  performed.  In  the  quarter  symmetric  3Z)  models  (Figure  9),  X-  and  Z-axes 
are  along  the  peened  surface  and  Y-axis  is  represent  depth  of  the  component 
(in  the  direction  of  the  UP  pressure).  In  the  biased  mesh  study,  the  element  size 
along  X-  and  Z-axes  is  kept  constant.  The  element  size  is  biased  along  the  Y- 
axis.  In  this  direction,  the  finest  element  is  at  the  top  surface.  The  biased  mesh 
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Figure  12:  Mesh  convergence  results 
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Table  3:  Mesh  convergence  study 


Number  of  Elements 

Mesh 

CPU  Time 

X-axis 

Y-axis 

Z-axis 

Type 

in  minutes 

40 

40 

40 

Uniform 

60 

80 

50 

80 

Biased 

720 

80 

60 

80 

Biased 

1,260 

80 

70 

80 

Biased 

1,680 

80 

80 

80 

Uniform 

1,950 

160 

160 

160 

Uniform 

10,800 

is  decided  using  two  parameters:  the  size  of  the  first  element  on  the  surface 
and  the  total  number  of  elements.  In  this  study,  the  size  of  the  top  element  is 
fixed  at  0.0625  mm,  which  is  equivalent  to  having  80  elements  per  5  mm  along 
Y-axis.  The  number  of  elements  used  is  50,  60,  or  70.  The  mesh  in  the  other 
two  directions  is  uniform  (80  elements). 

The  time  taken  by  each  of  these  analyses  is  shown  in  Table  3.  The  CPU  time 
shown  in  the  table  is  in  minutes  on  the  Glenn  cluster  at  OSC.  This  time  does 
not  include  the  waiting  time.  Typical  waiting  time  for  40  x  40  x  40  and  160  x 
160  X  160  meshes  is  around  15  and  2000  minutes,  respectively.  This  time  is 
not  the  exact  time.  The  time  in  minutes  is  rounded  to  the  nearest  10. 

The  extensive  mesh  convergence  study  provides  the  following  results,  indicat¬ 
ing  that  it  is  possible  to  achieve  convergence  in  3D  but  it  is  not  computationally 
feasible.  Comparison  of  the  results  shows  that  the  biased  mesh  is  able  to  pro¬ 
duce  similar  results  without  having  to  use  a  finer  mesh  for  the  entire  depth. 
A  biased  mesh  of  80  x  60  x  80  would  be  an  acceptable  mesh  for  research 
purposes.  The  above  convergence  study  is  performed  on  the  residual  stress 
profile.  In  the  design  optimization  of  the  LP  process  four  performance  metrics 
are  used.  These  metrics  are  calculated  form  the  residual  stress  profile.  The 
objective  of  all  optimization  formulations  is  to  compressive  stress  volume.  A 
convergence  study  is  also  be  performed  with  respect  to  the  compressive  vol¬ 
ume.  The  results  of  the  study  are  shown  in  table  4.  The  table  shows  that  com¬ 
pressive  volume  values  are  comparable  to  each  other.  This  study  shows  that 
the  selected  biased  mesh  can  be  assumed  as  converged  mesh  for  the  research 
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Table  4:  Mesh  Convergence  Study  with  Respect  Compressive  Volume 


Number  of  Elements 

Compressive 

X-axis 

Y-axis 

Z-axis 

Volume  (mm^) 

40 

40 

40 

4.17 

80 

50 

80 

4.48 

80 

60 

80 

4.51 

80 

70 

80 

4.53 

80 

80 

80 

4.63 

160 

160 

160 

4.67 

purposes. 

2.4.2  Preliminary  results 

The  LP  pressure  pulse  inputs  energy  into  the  material,  and  the  energy  is  di¬ 
vided  into  different  types  of  energies.  Part  of  the  energy  is  used  in  plastic 
deformation,  and  part  of  it  is  converted  into  kinetic,  internal,  and  heat  ener¬ 
gies.  Figure  13  shows  the  history  of  the  total  energy,  the  kinetic  energy,  the 
internal  energy,  and  the  plastic  dissipation  energy  for  a  3  GPa  peak  pressure 
pulse.  The  time  duration  of  2500  ns  for  dynamic  analysis  is  selected  based  on 
the  variation  in  the  plastic  deformation  and  the  internal  and  kinetic  energies. 
When  these  three  quantities  become  approximately  stable  at  2500  ns,  the  next 
step  of  the  analysis  begins. 

The  explicit,  implicit,  and  total  analysis  times  for  each  model  are  different. 
These  depend  upon  the  geometry,  number  of  elements,  analysis  duration,  and 
type  of  elements.  For  example,  a  model  with  infinite  elements  requires  less 
time  compared  to  a  model  without  infinite  elements.  The  reason  behind  this 
is  that  infinite  elements/region  tend  to  absorb  the  peening  energy  and  allows 
the  model  to  reach  equilibrium  state  faster.  In  the  absence  of  infinite  elements 
the  energy  reflects  back  from  the  geometric  boundaries.  In  such  cases,  the 
equilibrium  is  achieved  more  slowly.  Figure  13  shows  that  the  plastic  defor¬ 
mation  lasts  for  only  a  small  duration  when  the  applied  pressure  is  above  the 
elastic  limit.  The  figure  shows  that  the  external  work  and  the  plastic  dissipa¬ 
tion  energy  change  occur  for  the  duration  of  the  pressure  pulse,  which  is  170 
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Figure  13:  History  of  internal  energy,  artificial  strain  energy,  kinetic  energy,  plastic  dissipation  and  external  work 
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Figure  14:  Resultant  residual  stress  distribution  plot  for  3  GPa  peak  pressure  pulse 

ns.  After  that,  the  plastic  dissipation  energy  is  constant  at  approximately  4.5 
mJ.  The  kinetic  energy  increases  to  13.5  mJ  pressure  is  applied,  then  drops 
to  6.5  mJ  at  1000  ns,  and  goes  to  almost  zero  after  1500  ns.  There  are  many 
activities  going  on  during  the  initial  1000  ns,  and  all  energies  reach  an  approx¬ 
imately  stable  level  from  1000  ns  to  2000  ns.  It  can  be  inferred  that  there  is  no 
interaction  between  various  stress  waves  inside  the  component  after  2500  ns. 


Figure  14  shows  the  residual  stress  distribution  of  a  3  GPa  pulse  of  duration 
170  ns.  The  maximum  compressive  stresses  are  at  the  surface.  On  the  surface, 
the  location  of  the  maximum  stress  is  not  at  the  center  of  the  spot.  For  3  GPa 
location  of  maximum  compressive  stress  is  0.05  away  from  the  center.  The 
maximum  compressive  stress  location  is  not  independent  of  the  LP  parame¬ 
ters.  For  example  the  location  changes  with  the  pressure  pulse  magnitude. 
For  higher  pressure  magnitudes  the  location  moves  away  from  the  center  and 
below  the  top  surface.  The  volume  of  the  compressive  stress  is  less  as  com¬ 
pare  to  the  volume  of  tensile  stress.  The  maximum  tensile  stress  are  below 
the  compressive  stress  at  the  center  of  the  model.  The  depth  of  the  compres¬ 
sive  stress  is  approximately  same  in  the  entire  spot.  The  maximum  depth  of 
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tensile  stress  is  at  the  center  of  the  spot.  This  figure  shows  that  LP  is  able 
to  generate  compressive  stress  in  the  component,  but  there  are  two  causes  for 
concern.  With  the  current  set  of  parameters,  a  part  of  the  component  has  ten¬ 
sile  stresses.  Tensile  stresses  are  generated  to  balance  the  compressive  stress 
for  equilibrium  in  the  absence  of  external  force.  These  tensile  stresses  cannot 
be  eliminated;  however,  these  can  be  constrained  in  optimization  to  avoid  ad¬ 
verse  consequences.  Another  cause  for  concern  is  a  small,  shallow,  reduced 
compressive  or  tensile  stress  region  in  the  center  of  the  laser  spot  for  higher 
peak  pressures;  however,  it  is  not  representative  of  most  of  the  laser  spot. 

2.4.3  Experimental  Validation 

In  this  section,  finite  element  results  are  compared  with  the  experimental  re¬ 
sults  from  the  literature.  Multiple  analyses  are  performed  for  this  purpose.  The 
computed  residual  stresses  are  compared  with  the  residual  stress  field  from  the 
available  experimental  results.  Two  types  of  residual  stresses  are  reported: 
maximum  principal  and  uni-directional  stress  (Gh  or  G33).  As  in  the  LP  liter¬ 
ature,  most  calculations  in  the  report  are  use  the  uni-directional  stresses. 


Residual  Stress  Averaging 

All  the  experimental  residual  stresses  reported  here  are  taken  from  a  PhD  dis¬ 
sertation  by  Nam  [91  at  the  Ohio  State  University.  These  experimental  results 
are  obtained  using  the  X-ray  diffraction  technique.  This  is  a  destmctive  resid¬ 
ual  stress  measurement  technique  in  which  a  thin  layer  of  material  is  removed 
using  electro-polishing.  The  residual  stresses  obtained  using  X-ray  diffraction 
are  averaged  over  a  rectangular  area.  This  rectangular  area’s  size  and  shape 
are  dependent  on  the  size  and  shape  of  the  X-ray  diffraction  spot.  For  a  direct 
comparison  between  FE  simulation  results  of  residual  stress  distributions,  the 
same  averaging  should  be  used.  For  simulation  results,  the  averaging  is  done 
for  1  mm  square.  For  a  node,  residual  stress  values  at  all  the  nodes  in  0.5  mm 
distance  in  all  four  direction  are  added.  This  value  is  divided  by  the  total  num¬ 
ber  of  nodes.  The  resulting  value  is  the  residual  stress  value  for  the  selected 
node.  Figure  15  can  be  used  to  understand  the  process  averaging  process.  The 
figure  shows  a  layer  of  nodes  in  a  model.  The  red  node  shown  in  the  figure 
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1  mm 


Node 


Figure  15:  Averaging  Schematic 


is  selected  to  show  the  averaging  process.  The  residual  stress  values  of  the 
nodes  in  the  1  mm  square  are  added.  This  value  is  divided  by  total  number  of 
nodes  (25  in  this  figure).  The  resulting  value  is  the  residual  stress  value  at  the 
selected  red  node.  Similar  can  be  used  to  calculate  a  averaged  residual  stress 
profile  for  a  entire  model. 

This  is  necessary  to  add  that  the  same  averaging  process  is  used  for  all  the 
calculation  in  this  research.  This  process  is  independent  of  the  LP  parame¬ 
ters.  If  there  are  high  stress  concentration  location  (compressive  or  tensile), 
these  values  will  reduce  or  may  even  disappear  due  to  averaging.  For  a  re¬ 
search  that  is  interested  in  the  local  effects  of  LP,  averaging  may  not  be  a  good 
option.  Averaging  is  a  necessary  for  comparison  between  experimental  (with 
XRD  measurements)  and  simulation  results.  The  same  averaging  process  is 
used  in  this  research  for  consistency  purposes.  Depending  upon  the  research 
requirements  this  averaging  process  can  be  changed. 
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Result  Comparison 


In  each  case,  residual  stress,  strain,  and  deformation  properties  of  the  control 
volume  are  transfered  to  the  next  phase,  and  the  infinite  element  section  is 
remodeled.  The  symmetric  3D  model  is  analyzed  for  5.5  GPa  and  8.3  GPa 
pressure  pulses.  Results  obtained  for  5.5  GPa  and  8.3  GPa  are  compared  with 
published  experimental  [91  results.  The  comparison  of  the  simulation  and  ex- 


Figure  16:  Comparison  of  experimental  results  with  simulation  at  1  mm  from  spot  center  line 

perimental  results  shows  agreement  between  the  two  in  majority  of  the  com¬ 
ponent  depth.  The  simulations  are  able  to  estimate  the  depth  of  compressive 
stress  and  major  portion  of  compressive  stress  magnitude.  The  simulation  esti¬ 
mate  of  residual  stress  at  the  surface  show  difference  in  magnitude.  The  differ¬ 
ence  in  the  results  can  be  due  to  the  material  model,  variation  in  pressure  pulse 
properties,  and  uncertainties  in  the  experimental  measurements.  Although  im¬ 
provement  is  necessary  and  possible,  it  may  be  difficult  to  exactly  match  the 
experimental  results. 
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Figure  17:  Schematic  of  sequential  LP  at  multiple  locations  on  a  FEA  model 


2.4.4  Residual  Stress  Profile 

Figure  16  compares  the  experimental  and  simulation  results  for  the  average 
residual  stresses  of  the  two  peak  pressure  pulses,  5.5  GPa  and  8.3  GPa,  for 
one-sided  single  LP.  These  comparisons  are  made  along  a  line  offset  1  mm 
from  the  centerline  of  the  laser  spot.  Comparisons  are  made  at  this  location 
because  the  highly  localized  effect  of  the  release  waves  focuses  at  the  center 
of  the  spot,  causing  reduction  in  the  compressive  stress. 

2.5  Two  Shots  Sequence  Model 

To  investigate  the  effect  of  sequential  peening  and  overlap,  a  simple  exten¬ 
sion  of  quarter  symmetric  3Z)  model  is  developed.  As  shown  in  Figure  17, 
this  model  is  half  symmetric.  The  finite  element  properties  and  the  LP  prop- 


50 


erties  for  this  model  are  same  as  that  of  the  quarter  symmetric  3Z)  model.  The 
simulation  of  sequential  LP  at  multiple  locations  is  performed  for  two  config¬ 
urations.  In  the  first  configuration,  two  LP  locations  that  contact  each  other 
without  overlap  are  considered.  In  the  second  configuration,  there  is  a  50% 
overlap  between  the  two  LP  locations.  Figure  17  shows  a  schematic  of  se¬ 
quential  LP  at  multiple  locations  with  50%  overlap  between  the  two  locations. 
In  the  figure,  the  control  volume  has  two  overlaping  spots  and  infinite  elements 
all  around.  In  the  case  of  0%  overlap,  the  control  volume  dimensions  are  5  mm 
X  5  mm  X  15  mm,  and  in  the  case  of  50%  overlap,  the  dimensions  are  5  mm 
X  5  mm  X  12.5  mm.  Different  sizes  of  control  volumes  are  used  to  minimize 
computational  expenses. 

The  computational  time  for  analysis  is  huge.  Using  a  Pentium(R)  4  CPU, 
3.40  GHz  personal  computer  with  2.00  GB  of  RAM,  it  takes  approximately 
48  Hours  of  CPU  time  to  finish  one  analysis.  Typically,  depending  upon  the 
size  and  requirements  of  a  component,  it  may  need  from  20  to  100  LP  shots. 
Therefore,  the  time  taken  for  the  above  case  is  an  indication  that  FE  simulation 
for  a  practical  component  will  be  a  computational  challenge. 

2.6  Seven  Shots  Sequence  at  Multiple  Locations  Model 

In  this  section,  a  cylindrical  model  is  developed  that  can  be  employed  to  inves¬ 
tigate  seven  shot  sequence.  The  model  is  shown  in  Figure  18.  The  cylindrical 
shape  is  selected  to  minimize  the  computational  requirements  of  a  seven  shot 
sequence.  The  cylindrical  shape  is  better  suited  to  the  investigation  compared 
to  a  rectangular  component.  This  model  requires  a  lower  number  of  elements 
compared  to  a  rectangular  model.  The  control  volume  of  the  model  has  finer 
mesh.  The  rest  of  model  has  coarser  mesh.  The  model  can  also  be  employed 
to  investigate  the  effect  of  multiple  shots  in  partial  overlapping  conditions. 

2.7  Parametric  Plate  Model 

The  parameters  of  validated  symmetric  3D  simulation  methodology  are  used 
to  develop  the  parametric  plate  model.  Figure  19  shows  the  FE  model  of  the 
plate  with  a  control  volume  size  oflOx  10x5  (mrn^).  The  model  is  a  infinite 
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Figure  18:  Cylindrical  model  for  sequence  investigation 
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Figure  19:  FEA  model  of  a  plate  with  five  shots 
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slab  because  of  the  presence  of  infinite  elements.  In  this  report,  it  is  called  a 
parametric  plate  because  it  can  consider  most  of  the  LP  parameters.  As  shown 
in  the  figure,  infinite  elements  extend  from  the  plate  by  5  mm  on  four  sides 
and  on  the  bottom  of  the  plate.  In  this  model,  up  to  five  different  locations  can 
be  selected  for  peening.  The  figure  shows  four  spot  locations;  a  fifth  one  is 
located  at  the  center  of  the  plate.  These  locations  are  variables  to  generate  the 
desired  overlapping  configurations  or  totally  separate  spots.  Apart  from  the 
parametric  locations  (or  overlap),  the  model  allows  the  designer  to  specify  the 
spot  size  and  shape,  pressure  pulse  magnitudes,  shapes  and  duration,  number 
of  shots  at  the  same  location,  and  the  sequence  of  shots.  This  model  provides 
many  advantages,  but  introduces  the  issue  of  computational  resources.  For  the 
above-mentioned  methodology,  a  five  shot  simulation  (at  the  same  location  or 
different)  takes  4  days  on  the  Glenn  cluster  at  OSC.  Design  optimization  is  an 
iterative  process  and  the  use  of  metamodels  is  advised  to  avoid  prohibitively 
large  computational  expense. 

2.8  Rectangular  Coupon  Model 

A  rectangular  coupon  is  selected  for  modeling  because  it  can  be  employed  in 
future  work  for  experimental  validation  of  residual  stress  profiles  and  fatigue 
life  estimation.  The  basic  framework,  mesh  convergence,  and  material  model 
for  the  coupon  are  taken  from  the  validated  3D  model.  The  region  of  interest  is 
called  the  control  region.  The  controlled  region  is  uses  a  very  fine  mesh.  The 
region  surrounding  the  control  region  is  meshed  using  coarse  mesh.  These 
region  can  be  clearly  seen  in  Figure  20.  The  coupon  model  can  be  used  to 
perform  laser  peening  simulations  to  estimate  induced  residual  stresses.  The 
same  coupon  with  residual  stresses  can  be  employed  to  apply  fatigue  loading. 
This  setup  will  allow  the  calculation  of  the  fatigue  life  of  the  peened  coupon, 
and  the  ability  to  compare  it  with  the  fatigue  life  of  the  unpeened  coupon.  For 
the  coupon  model,  a  three  shot  simulation  takes  8  CPU  hours  on  a  processor 
at  OSC. 
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-  Test  Loading 


Figure  20:  FE  model  of  the  four  point  bend  specimen 


2.9  Chapter  Summary 

This  section  began  with  an  explanation  of  the  simulation-based  methodology. 
The  methodology  is  adopted  to  predict  residual  stresses  generated  by  the  LP 
pressure  pulse.  Based  on  the  literature,  a  range  of  LP  parameters  is  determined 
and  used  for  simulation.  Initially,  this  methodology  is  successfully  verified  by 
comparing  the  simulation  results  with  the  available  experimental  results.  Then 
the  methodology  and  model  are  used  to  develop  models  to  investigate  two 
overlapping  shots,  a  seven  shot  sequence,  and  parametric  design  optimization 
of  a  laser  peening  process.  In  summary,  this  work  added  to  the  understand¬ 
ing  of  the  physical  process  of  laser  peening  and  developed  a  simulation-based 
methodology,  which  can  be  used  for  analysis  of  complex  geometries  repre¬ 
senting  practical  applications.  Parametric  investigations,  residual  stress  field 
optimization,  and  residual  stress  estimation  for  fatigue  life  calculation  can  be 
performed  by  using  the  simulation  methodology  and  the  finite  element  models 
presented  in  this  section. 
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3.0  LP  OPTIMIZATION:  INTRODUCTION  AND  PARAMETRIC 

INVESTIGATIONS 

Mathematical  optimization  is  a  process  of  identifying  the  best  suitable  design 
from  a  collection  of  alternatives  without  having  to  experiment  all  possible  al¬ 
ternatives.  There  are  many  simple  problems  in  which  optimal  parameters  can 
be  obtained  through  experience  and  experiments.  However,  for  a  nonlinear 
time-dependent  mixed- variable  process  such  as  LP,  it  is  difficult  determine  op¬ 
timal  parameters  using  this  approach.  It  is  necessary  to  perform  a  systematic 
optimization  of  the  LP  process. 

The  LP  process  has  many  variables,  such  as  the  pressure  pulse  magnitude  and 
mid-span,  spot  shape  and  size,  the  location  and  layout  of  the  spots,  and  the 
number  of  shots  at  the  same  location.  For  a  comprehensive  optimization, 
parametric  investigations  are  needed  to  determine  the  performance  function 
sensitivities  with  respect  to  optimization  variables.  The  optimization  of  the 
LP  process  is  a  mixed- variable  optimization  problem  because  pressure  pulse 
magnitude  and  duration,  spot  shape  and  size,  and  amount  of  overlap  are  con¬ 
tinuous,  while  the  number  of  shots  and  spot  layout  are  discrete.  The  time 
dependent  nature  of  the  process  increases  the  complexity  because  each  shot  is 
a  nonlinear  plastic  process.  This  process  creates  residual  stresses  and  proper¬ 
ties  of  these  stresses  affects  the  selection  of  the  parameters  of  the  subsequent 
shots. 

Overall,  the  main  objectives  of  this  section  are  (1)  to  present  the  mathemat¬ 
ical  formulation  of  a  mixed- variable  optimization  problem,  (2)  to  introduce 
gradient-  and  non-gradient-based  optimization  methods  and  surrogate  models, 
and  (3)  to  perform  parametric  investigations. 

3.1  Mathematical  Eormulation  of  Mixed- variable  Optimization 

Sometimes,  an  optimization  problem  has  input  variables  that  are  continuous, 
integer,  and  discrete.  Such  problems  are  called  mixed-variable  optimization 
(MVO)  problems.  A  typical  mathematical  formulation  of  an  MVO  problem  is 
defined  as 
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Maximize:  /(^) 

Sub ject  to  : 

Equality  constraints  :  hi{X)  =0,  i=  1,2, 
Inequality  constraints gi{X)<Q,  /  =  m+ l,m  +  2, 

The  input  vector  is 
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••5 

\ 
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The  variable  bounds  on  each  input  variable  are  given  by 


/=  l,2,...,nc, 

xf  <  x\  <  xf  ,  i  =  1,2, 
xf<xf<xf“,  i  =  l,2,...,nj. 

/  =  1,2, 

where  /(2C)  is  the  objective  function,  xf  G  7?^,  x-  E  R\xf  E  R‘^,  and  x^  E  R^ 
denote  sets  of  continuous,  integer,  discrete,  and  binary  (zero-one)  variables, 
respectively,  and  xf  and  x*“  denote  the  lower  and  upper  bounds  of  the  re¬ 
spective  variables  (replace  (*)  with  c,i,d,  orb).  The  total  number  of  vari¬ 
ables  is  n  =  Uc  +  ni  +  uj  +  nt,  and  the  total  number  of  constraints  is  p.  Here 
nc,ni,n^,  and  are  the  number  of  continuous,  integer,  discrete,  and  binary 
variables. 

3.2  Optimization  Methods 

Optimization  algorithms  are  broadly  divided  into  two  categories:  (i)  Gradient- 
based  algorithms,  and  (ii)  Non-gradient-based  algorithms.  All  the  optimiza¬ 
tion  algorithms  are  iterative  in  nature.  The  gradient-based  algorithms  move 
from  the  current  point  to  the  next  point  based  on  gradients/sensitivities.  In  this 
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process  of  moving  points,  the  objective  function  and  constraint  values  tend  to 
improve.  In  a  population-based  approach,  in  each  iteration,  based  on  a  set  of 
rules,  a  new  population  is  generated  from  the  previous  population.  Population 
based  methods  may  not  use  gradient  information. 

There  are  many  algorithms  (both  gradient  and  non-gradient)  that  treat  all  kinds 
of  variables  as  continuous  variables  and  round-off  the  values  of  integer  and/or 
discrete  variables.  Some  methods  also  use  the  penalty  function  approach  to 
handle  integer  and  discrete  variables.  Most  gradient-based  algorithms  are  de¬ 
signed  for  optimization  problems  with  only  continuous  variables  because  of 
their  need  for  gradients.  An  algorithms  to  solve  such  problems  is  branch 
and  bound  [741.  Recently,  many  evolutionary  techniques  (genetic  algorithm 
and  particle- swarm  optimization)  have  been  used  to  solve  mixed- variable  op¬ 
timization  problems. 

3.2.1  Gradient-Based  Methods 

Gradient-based  optimization  methods  use  the  objective  function,  constraints 
functions,  and  gradient  information  to  guide  the  search  strategy.  The  gradient 
information  can  be  obtained  from  the  first-  and/or  second-order  derivatives 
of  the  objective  function  and  constraints.  The  formulation  of  gradient-based 
optimization  is  given  by  Equation  20: 


#+i)  (20) 

where  is  the  next  point,  is  the  current  point,  is  the  step-length 

parameter  or  distance  to  travel  to  reach  the  (k+  point  from  and 
is  the  unit  direction  vector.  Most  gradient  based-optimization  algorithms  are 
about  finding  the  direction  vector  and  distance  (a*^^^)  to  travel  in  that 
direction. 

Popular  methods  to  find  the  appropriate  direction  are  Newton’s  method,  the 
Marquardt  modification  of  Newton’s  method,  the  steepest  descent  method, 
the  gradient  projection,  the  reduced  gradient  method,  the  conjugate  gradient 
method,  the  Davidon-Fletcher-Powell  method,  the  feasible  direction  method, 
and  the  penalty  function  methods.  Many  of  the  direction  finding  algorithms 
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mentioned  above  also  have  approaches  to  find  ,  and  other  popular  methods 
use  simple  differentiation,  the  equal  interval  search,  the  bisection  method,  and 
the  golden  section  search. 

3.2.2  Non-gradient  Based  Optimization 

Generally,  non-gradient  based  optimization  algorithms  are  population-based 
probabilistic  methods,  also  called  Evolutionary  Algorithms  (EA).  These  al¬ 
gorithms  may  not  use  gradient  information.  Simulated  Annealing  (SA),  Ge¬ 
netic  Algorithms  (GA)  [75,761,  Ant  Colony  Optimization  (AGO),  and  Particle 
Swarm  Optimization  (PSO)  are  a  few  examples.  Two  of  these  methods,  the 
Genetic  Algorithm  and  Particle  Swarm  Optimization  are  briefly  presented  in 
the  following  sections. 


Genetic  Algorithm 

GA  belongs  to  the  class  of  natural  adaptive  optimization  algorithms  that  mimic 
the  process  of  natural  evolution.  Since  the  pioneering  work  of  John  Holland, 
who  proposed  the  possibility  of  a  genetics-based  search  to  solve  optimiza¬ 
tion  problems,  genetic  algorithms  have  emerged  as  a  robust,  stochastic  search 
method  in  complex  problem  domains  ranging  from  process  optimization,  con¬ 
trols,  scheduling,  motion  planning,  pattern  recognition,  and  stmctural  opti¬ 
mization.  However,  since  GA  is  a  population-based  approach,  it  requires  a 
large  number  of  function  evaluations. 


Particle  Swarm  Optimization 

PSO  is  an  algorithm  proposed  by  Kennedy  and  Eberhart  [771,  motivated  by 
social  behavior  of  organisms  such  as  bird  flocking  and  fish  schooling.  Like 
evolutionary  algorithms,  PSO  is  a  population  based  approach.  This  population 
is  modified  from  the  current  iteration  to  the  next  iteration  based  on  a  set  of 
rules.  This  set  of  rules  tends  to  mimic  the  problem  solving,  social  behavior  of 
group  of  birds.  In  this  set  of  rules  each  individual  (particle)  has  its  fitness,  and 
a  communication  structure,  assigning  neighbors  for  each  individual  to  interact 


58 


with.  These  individuals  iteratively  calculate  their  fitness  solutions  and  remem¬ 
ber  the  location  of  their  best  fitness  so  far.  The  information  of  the  best  locations 
of  each  neighboring  individuals  is  available  to  all  neighboring  individuals  and 
the  best  of  all  the  individuals  is  also  available  to  every  individual.  Based  on 
this  information,  movement  parameters  of  each  individual  particle  are  evalu¬ 
ated.  These  movement  parameters  change  over  the  course  of  iterations  and  the 
particles  tends  to  converge  and  provide  an  optimum  solution.  Although  PSO 
is  a  powerful  algorithm,  it  tends  to  be  converge  slowly  on  some  problems,  and 
as  a  result  takes  a  large  number  of  function  evaluations. 

3.2.3  Advantages  and  Disadvantages 

This  section  discusses  the  advantages  and  disadvantages  of  gradient-  and  non¬ 
gradient-based  methods.  The  gradient-based  methods  tend  to  need  a  large 
number  of  function  evaluations  for  gradient  evaluation.  EAs  tend  to  be  robust, 
but  also  require  a  large  number  of  function  evaluations.  This  is  a  major  issue 
for  problems  like  LP  because  of  the  humongous  computational  time  and  space 
required  to  perform  a  large  number  of  simulations.  The  disadvantage  of  gra¬ 
dient  methods  is  that  the  final  solution  may  depend  on  the  starting  point  of  the 
algorithm  and  tends  to  converge  on  local  optima.  Nevertheless,  these  methods 
clearly  improve  next  point  and  termination  for  an  optimal  solution.  EAs  have 
a  disadvantage:  a  solution  in  the  population  is  always  compared  to  other  solu¬ 
tions  in  the  selection,  not  in  absolute  terms,  and  these  methods  do  not  have  the 
concept  of  optimal  solution. 

EAs  can  approach  a  near  optimal  solution  quickly,  but  can  take  a  large  num¬ 
ber  of  iterations  to  find  an  exact  optimum.  EAs  are  suitable  for  parallel  pro¬ 
cessing  environment  because  of  the  population-based  approach,  but  gradient 
approaches  are  not  because  the  calculation  of  the  next  point  depends  on  the 
previous  point  (Equation  20).  However,  parallel  computing  can  be  utilized 
to  calculate  derivatives/gradients  for  gradient-based  methods.  Gradient-based 
methods  are  not  suitable  for  handling  discrete  variables  without  adjustment, 
but  evolutionary  algorithms  can  easily  be  adjusted  to  handle  discrete  variables. 

PSO  is  selected  to  solve  the  LP  optimization  problems.  Along  with  all  the 
advantages  of  EAs,  PSO  is  selected  because  of  its  ease  of  implementation. 
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Figure  21:  Flow-chart  of  the  simulation-based  optimization  process 


lower  number  of  user  parameters,  and  the  use  of  historical  information. 

3.3  Automated  Data  Transfer  Procedure 

Due  to  computing  and  licensing  limitations,  the  computing  is  done  in  two 
parts.  The  simulation-based  design/  optimization  process  generates  the  input 
FEA  files  on  a  personal  computer  (Windows  PC).  The  LP  finite  element  analy¬ 
ses  are  performed  at  OSC  (Unix).  To  perform  analysis,  data  must  be  transfered 
between  these  to  computers.  A  manual  data  transfer  method  is  not  the  most 
efficient  option.  An  automated  procedure  is  developed  that  does  not  require 
login  and  password  for  each  transfer  from  PC  to  supercomputer.  Although  the 
process  does  not  require  password,  it  is  secure  as  per  current  security  stan¬ 
dards.  A  one  time  paraphrase  is  used  for  security.  A  modern  programming 
language,  PuTTy,  is  used  for  automated  data  transfer.  A  successful  transfer 
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depends  upon  the  reliability  of  the  connecting  network  between  computers. 
However,  the  procedure  is  independent  of  minor  network  issues.  This  inde¬ 
pendence  is  achieved  by  iterative  attempts.  Many  attempts  are  made  to  transfer 
data.  The  time  between  two  attempts  increases  exponentially  with  the  number 
of  failed  attempts.  An  optimization  loop  will  not  crash  if  a  PC  is  not  able  to 
connect  to  the  server  for  few  minutes  or  few  attempts. 

Access  to  OSC,  knowledge  of  basic  Unix  and  batch  processing,  and  experi¬ 
ence  in  a  programming  language  are  a  few  requirements  that  must  be  met  to 
successfully  implement  the  procedure.  The  software  packages  MATLAB  and 
PuTTY  are  employed  to  develop  the  procedure  (Figure  21).  The  activities  on 
the  left  hand  side  of  the  vertical  line  in  the  figure  are  performed  on  a  PC.  The 
activities  on  the  right  hand  side  of  the  line  are  performed  at  the  OSC.  These 
software  packages  provide  a  set-up  that  generates  input  files  on  a  PC,  trans¬ 
fers  the  files  to  OSC,  performs  the  analyses  at  OSC,  and  brings  result  files  to 
the  PC.  A  brief  desorption  of  the  procedure  is  given  in  the  following  sections. 
Further  details  can  be  seen  in  Ref.  [781. 

3.3.1  One-time  Setup 

There  is  a  one-time  setup  needed  to  use  the  process.  The  setup  includes  two 
steps,  installing  PuTTy  and  setting  up  a  password-less  SSH  connection.  In  the 
first  step,  PuTTy  software  is  installed.  The  software  is  a  legally  free  to  down¬ 
load.  This  software  is  to  implement  Telent  and  SSH  for  Win32  and  Unix  plat¬ 
form.  In  the  second  step,  a  password  less  paraphrase  protected  process  is  in¬ 
stalled.  The  details  of  the  process  installation  are  available  at  the  OSC  website 
(http://www.osc.edu/supercomputing/software/apps/mdcs.shtml).  Only  the  step 

2. 1  is  needed  from  the  web-page.  The  remaining  steps  on  the  web-page  are  for 
Parallel  Computing  Toolbox  (PCT)  of  MATLAB  software.  For  assistant  and 
any  issue,  OSC  user  can  contact  OSC  help  at  oschelp@osc.edu.  After  one¬ 
time  setup  is  complete,  the  procedure  is  verified  as  given  by  at  the  end  of  the 
same  step  2. 1 . 
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3.3.2  Procedure  Steps 


In  this  research,  the  main  organizing  tool  of  the  procedure  is  MATLAB.  Any 
programming  language  such  as  C,  C++,  C#,  fortran,  python,  or  VB  can  be 
used  for  the  same  purpose.  The  optimization  code  (PSO),  finite  element  anal¬ 
ysis  input  files  generation  codes  (for  ABAQUS),  the  PuTTY  management,  and 
fatigue  analysis  files  generation  code  (for  fe-Safe)  are  written  in  MATLAB. 
As  shown  in  Figure  21,  the  procedure  involves  following  steps. 

(a)  Generate  finite  element  input  files  at  the  PC.  The  number  of  files  depend 
upon  the  number  of  particles  in  the  optimization. 

(b)  Employ  PuTTy  to  transfer  all  files  to  the  OSC  cluster. 

(c)  Perform  FEA  at  OSC  using  ABAQUS.  The  request  to  start  all  analyses  is 
given  at  the  same  time. 

(d)  Generate  the  required  result  files. 

(e)  Transfer  the  results  files  from  the  OSC  cluster  to  the  PC  using  PuTTy. 

(f)  Perform  additional  analysis  at  PC,  if  required. 

(g)  Check  for  the  convergence.  Stop,  if  converged.  If  not  converged,  create 
new  points  from  the  previous  points  and  the  results.  Repeat  the  process. 

3.3.3  Advantages  of  the  Procedure 

This  procedure  overcomes  the  licensing  limitations  and  optimally  utilizes  the 
computing  facility  of  OSC.  This  procedure  is  ideally  suited  for  the  research 
groups  that  do  not  have  supercomputers  available  at  hand  and  depend  upon 
OSC.  Overall,  the  automation  procedure  has  many  advantages: 

(i)  This  procedure  is  suitable  for  a  population-based  optimization  methods 
(PSO).  The  procedure  assists  in  performing  many  analyses  (equal  to  pop¬ 
ulation  size)  at  the  same  time  at  OSC. 

(ii)  Many  design  processes  involve  more  than  one  software  packages  to  com¬ 
pute  the  performance  function  (ABAQUS  and  fe-Safe  for  fatigue  life  op¬ 
timization  using  laser  peening).  The  developed  procedure  removes  the 
limitations  that  all  the  licenses  must  be  available  on  a  particular  PC/cluster. 
This  is  due  to  the  automated  data  transfer  between  computers. 
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(iii)  This  procedure  allows  an  optimal  utilization  of  the  OSC  computing  and 
licensing  resources. 

3.4  Surrogate  Models 

A  surrogate  model  or  function  approximation  is  a  mathematical  model  that 
represents  the  relations  between  dependent  variables  or  responses  and  inde¬ 
pendent  variables.  Surrogate  models  play  an  important  role  in  the  design  of 
experiments  as  models  for  prediction,  in  process  optimization,  and  in  un¬ 
certainty  quantification,  especially  for  large-scale  analysis.  There  are  many 
techniques  of  approximation,  such  as  one-point,  two-point,  and  multi-point 
models.  Among  one-point  surrogate  models,  linear,  reciprocal,  and  conser¬ 
vative  metamodels  are  the  most  commonly  used  in  structural  optimization. 
Two-point  Adaptive  Nonlinear  Approximation  (TANA),  Improved  Two-point 
Adaptive  Nonlinear  Approximation  (TANA2),  and  the  Response  Surface  Method 
(RSM)  are  a  few  examples  of  two-point  and  multi-point  metamodels. 

The  use  of  surrogate  models  in  general  as  well  as  in  this  research  is  not  to 
replace  the  full  scale  simulation  with  approximations.  Surrogate  models  are 
employed  during  intermediate  steps  in  an  optimization  process  with  the  aim 
of  reducing  the  computational  time.  The  full-scale  simulations  are  performed 
after  each  iteration.  The  goal  is  to  employ  an  effective  surrogate  model,  and 
a  minor  difference  in  a  surrogate  model  prediction  tends  not  to  affect  the  final 
solution.  The  effectiveness  of  the  surrogate  model  does  affect  the  number  of 
iterations  and  computational  time  required  to  reach  a  converged  solution  in  an 
iterative  algorithm. 

3.5  Parametric  Investigations 

The  work  discussed  in  section  3  provides  a  validated  simulation  procedure  for 
the  LP  process.  At  the  validation  stage  for  comparison  between  the  experi¬ 
mental  and  simulation  results,  the  parameter  selection  for  the  simulations  was 
constrained  by  the  available  experimental  results.  In  parametric  investigations, 
different  practically  possible  parameter  variations  are  explored  to  determine 
their  effect  on  the  residual  stress  field.  The  parametric  studies  employ  models 
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presented  in  section  3.  These  models  include  symmetric  3Z)  (Figure  9),  a  mod¬ 
ified  version  of  symmetric  3Z),  sequential  LP  simulation  model  (Figure  17), 
and  Cylindrical  model  (Figure  18).  There  are  multiple  objectives  and  multiple 
constraints  [791  in  the  LP  process.  Out  of  those,  four  performance  metrics- 
depth  of  compressive  stress  (CD),  maximum  compressive  (MC)  stress,  maxi¬ 
mum  tensile  stress  (MT),  and  total  compressive  volume  (CV)-are  selected  for 
comparison. 

These  performance  metrics  are  calculated  after  the  averaged  residual  stress 
profile  is  calculated  for  the  entire  model.  The  CD  is  calculated  at  a  selected  lo¬ 
cation.  The  number  of  nodes  (along  the  depth)  in  compression  at  the  selected 
location  are  determined.  The  number  of  elements  are  one  less  than  the  number 
of  nodes.  The  number  of  elements  multiplied  by  element  thickness  provides 
the  compressive  stress  depth  at  the  selected  location.  The  maximum  com¬ 
pression  and  tensile  stresses  are  determined  using  maximum  and  minimum 
finding  algorithm.  These  algorithms  find  the  values  from  the  entire  profile. 
The  compressive  volume  is  calculated  by  determining  the  number  of  nodes  in 
compression.  The  nodes  near  the  surface  are  checked  for  compressive  stresses. 
The  compressive  stress  region  that  is  generated  at  the  bottom  of  the  component 
(below  tensile  region)  is  not  considered  in  the  volume.  The  number  of  nodes  in 
compression  are  divided  by  the  total  number  of  nodes.  The  resulting  quantity 
is  then  multiplied  by  the  total  volume  in  compression. 

3.5.1  Temporal  Variation  of  Pressure 

This  parameter  is  defined  by  the  variation  of  pressure  with  time,  also  known  as 
pressure  pulse  shape.  The  pressure  pulse  shape  is  dependent  upon  the  tempo¬ 
ral  variations  of  the  laser  beam.  For  LP  applications,  two  types  of  laser  shapes 
are  typically  used.  The  first  type  of  shape  is  similar  to  a  Gaussian  distribution 
curve.  It  has  a  gradual  rise  and  similar  decay.  This  shape  takes  74.5  ns  to  reach 
its  peak  value,  stays  at  maximum  for  the  next  3  ns,  and  then  drops  to  zero  in 
the  next  74.5  ns.  In  the  second  shape,  the  pressure  increases  to  the  peak  within 
3  ns,  stays  at  maximum  for  the  next  3  ns,  and  then  reduces  slowly  to  zero  in 
152  ns.  The  areas  under  both  curves  are  the  same,  and  the  duration  is  152  ns. 
Both  shapes  are  shown  by  dashed  lines  in  Figure  22.  It  should  be  mentioned 
here  that  the  total  pulse  duration  is  152  ns  for  this  part  of  the  research;  other- 
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Figure  22:  Temporal  loading  profile  of  pressure  pulse 


wise,  pulse  duration  is  200  ns  (used  in  experimental  validation)  for  the  rest  of 
the  report.  The  peak  magnitude  in  this  study  is  5.5  GPa.  The  total  duration 
can  be  decreased  or  increased  depending  upon  the  laser  generation  mecha¬ 
nism.  The  sharp  and  slow  shapes  generate  compressive  depth  (mm)  of  1.46 
and  1.21,  maximum  compression  {MPa)  of  —478.6  and  —489.5,  maximum 
tension  {MPa)  of  139.0  and  127.8,  and  compressive  volume  (mm^)  of  7.14 
and  6.84,  respectively.  These  results  show  that  both  generate  approximately 
similar  results. 

The  magnitude  of  the  peak  pressure  depends  upon  the  laser  beam  intensity.  In 
addition,  the  use  of  a  confining  mechanism  increases  the  peak  pressure  by  4  to 
10  times.  Increasing  the  pressure  by  increasing  the  laser  beam  intensity  is  ul¬ 
timately  limited  by  the  onset  of  dielectric  breakdown  in  the  confinement  water 
layer  at  higher  laser  beam  intensities.  When  this  breakdown  occurs,  the  laser 
energy  no  longer  couples  into  the  expanding  plasma  underlying  the  water,  and 
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Table  5;  Results  of  pressure  magnitude  and  spot  radius  comparison 


Pressure  Magnitude  {GPa) 

Spot  Radius  ( 

mm) 

4.0 

6.0 

8.3 

1.00 

1.75 

2.50 

CD  (mm) 

0.61 

0.98 

1.6 

0.86 

0.97 

0.96 

CV  (mm) 

3.38 

5.87 

8.12 

1.42 

3.34 

5.17 

MC  {MPa) 

-282.2 

-487.3 

-831.3 

-212.4 

-362.2 

-374.6 

MT  {MPa) 

39.8 

79.94 

199.1 

36.04 

46.2 

60.4 

the  pressure  increase  ceases.  Although  simulations  show  that  both  temporal 
profiles  generate  similar  results,  to  avoid  the  effects  of  the  breakdown  phe¬ 
nomenon  a  sharp-rise  pulse  is  preferred  over  a  Gaussian  pulse.  In  the  case  of  a 
sharp  rise,  peak  pressure  is  built  in  the  very  beginning,  before  the  breakdown 
happens. 


3.5.2  Pressure  Pulse  Magnitude 

The  magnitude  of  the  pressure  pulse,  which  is  dependent  on  laser  beam  inten¬ 
sity,  is  the  most  important  parameter  in  the  LP  process.  Three  peak  pressure 
pulses  of  4.0  GPa,  6.0  GPa,  and  8.3  GPa,  are  used  to  compare  the  effects. 
The  pressure  pulse  shape  and  duration  are  shown  using  a  solid  line  in  Figure 
22.  The  results  are  shown  in  Table  5.  The  results  show  that  higher  pressure 
produces  higher  depth,  magnitude,  and  volume  of  compressive  stress,  and  an 
increase  in  the  tensile  magnitude. 

3.5.3  Spot  Size 

As  mentioned  in  Section  3.1.1,  a  bigger  spot  size  tends  to  produce  a  higher 
magnitude  and  volume  of  compressive  residual  stress  for  the  same  peak  pres¬ 
sure.  This  fact  is  verified  in  this  parametric  research,  in  which  three  spot  radii- 
1.0,  1.75,  and  2.5  mm-are  used  for  comparison  without  changing  any  other 
parameters.  The  pressure  pulse  magnitude  is  5.0  GPa.  The  four  performance 
metrics  are  shown  in  Table  5.  The  simulation  results  show  that  the  bigger  spot 
size  produces  higher  magnitude,  depth,  and  volume  of  compressive  stress.  In 
these  comparisons,  location  is  not  considered.  The  value  of  maximum  com- 
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Three  spot  shapes 


Figure  23:  Three  spot  shapes  investigated 


pression  can  be  different  in  the  three  cases.  The  increase  in  volume  can  be 
attributed  to  the  increase  in  the  size  of  the  spot  shape. 

3.5.4  Spot  Shape 

As  noted  in  Section  3.1.1,  rectangular  and  circular  shapes  are  used  most  often. 
This  section  investigates  only  small  variations  of  circular  shapes.  These  vari¬ 
ations  are  generated  by  changing  the  ratio  of  the  small  diameter  to  the  large 
diameter  of  an  ellipse.  The  three  spot  shapes  used  in  this  study  are  shown 
in  Figure  23.  This  figure  shows  a  quarter  of  each  spot  shape.  The  three  spot 
shapes  are  generated  by  considering  three  diameter  ratios  (1.0,  0.75,  and  0.50). 
The  pressure  pulse  magnitude  is  5.0  GPa.  The  radius  of  the  circle  (diameter 
ratio  1.0)  is  2.50  mm.  The  highest  diameter  of  remaining  elliptical  shapes 
(diameter  ratio  0.75  and  0.50)  is  2.50  mm.  The  smallest  diameters  are  1.875 
mm  (for  diameter  ratio  0.75,  0.75  *  2.50  =  1.875)  and  1.25  mm  (for  diame¬ 
ter  ratio  0.50,  0.50  *  2.50  =  1.25).  As  shown  by  the  solid  lines  in  the  figure 
for  the  diameter  ratio  1.0,  the  spot  takes  a  circular  shape.  The  residual  stress 
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Figure  24;  Residual  stress  field  for  three  spot  shapes 


fields  at  the  center,  shown  in  Figure  24,  are  compared  to  the  three  cases.  The 
figure  shows  a  reduction  in  the  compressive  stresses  at  the  surface,  with  in¬ 
crease  in  the  diameter  ratio.  For  the  ratios  1.0,  0.75,  and  0.50,  the  maximum 
compressive  stresses  are  —374.6,  —387.4,  and  —405.71,  respectively.  These 
represent  3.6%  and  8.5%  increase  from  the  baseline  values.  The  maximum 
tensile  stresses  are  60.4,  48.7,and  35.31,  respectively.  These  represent  19.4% 
and  41.5%  decrease  from  the  baseline  values.  These  results  indicate  that  spot 
shaoe  change  can  significantly  reduce  the  maximum  tensile  stresses  with  min¬ 
imal  change  in  the  maximum  compressive  stresses.  The  four  performance 
metrics  (depth  of  compressive  stress,  maximum  compressive  stress,  maximum 
tensile  stress,  and  total  compressive  volume)  for  each  case  are  shown  in  Table 
6.  The  table  shows  that  variation  in  compressive  and  tensile  magnitude  are  not 
similar.  The  reduction  in  tensile  stress  is  higher  as  compare  to  compressive 
stress.  This  parameter  should  be  employed  to  comtrol  the  tensile  stress. 
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Table  6;  Results  of  spot  shape  comparison 


Diameter  ratio 

1.00 

0.75 

0.50 

Depth  (mm) 

0.96 

1.1 

1.0 

Volume  (mm) 

5.17 

3.97 

2.56 

Maximum  Compression  {MPa) 

-374.6 

-387.4 

-405.71 

Maximum  Tension  {MPa) 

60.4 

48.7 

35.31 

3.5.5  Thickness  of  Component 

The  components  to  be  peened  come  in  various  shapes  and  sizes,  and  these 
boundary  conditions  can  significantly  affect  the  peening  process.  Therefore, 
component  thickness  is  a  parameter  that  requires  close  examination.  Since  the 
thickness  of  a  component  is  the  distance  traveled  by  the  shock  waves  in  one 
direction,  it  can  have  considerable  effect  on  shock  wave  propagations  and  the 
resulting  residual  stress  profile.  In  this  section,  the  effect  of  the  component 
thickness  is  investigated  for  cases  when  the  peening  surface  and  the  surface 
on  the  opposite  side  are  parallel  to  each  other.  The  FE  model  for  thickness 
investigations  is  shown  in  Figure  25.  The  model  is  a  modified  version  of  the 
quarter  symmetric  3Z)  model  shown  in  Figure  9.  The  model  is  modified  by 
removing  the  infinite  elements  in  the  bottom  of  the  model  and  replacing  them 
with  finite  elements.  The  total  thickness  of  the  component  varies  from  3  mm 
to  10  mm  in  increments  of  1  mm.  Normally  components  less  than  3  mm  in 
thickness  tend  not  to  have  parallel  surfaces.  The  element  size  changes  in  the 
direction  marked  by  thickness  in  Figure  25.  Using  the  thickness  and  mesh 
convergence  study,  an  appropriate  unbiased  mesh  is  selected.  A  finer  mesh  is 
used  for  smaller  thicknesses  (<  5  mm)  and  a  comparatively  coarser  mesh  is 
used  for  larger  thicknesses.  In  both  cases  the  mesh  is  equal  to  or  finer  than  the 
converged  mesh.  The  residual  stress  Gn  along  the  thickness  at  1  mm  from  the 
center  is  plotted  for  all  the  cases. 

A  comparison  of  residual  stress  profiles  for  various  thicknesses  is  shown  in 
Figure  26.  These  results  show  that  the  component  thickness  affects  the  resid¬ 
ual  stress  at  the  surface  when  the  thickness  is  less  than  4  mm.  The  residual 
stress  profile  in  the  area  of  interest  is  not  affected  significantly  by  the  thick- 
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Figure  25;  Residual  stress  field  comparison  for  depth  from  3  mm  to  10  mm 
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Depth 

Figure  26;  Residual  stress  field  comparison  for  depth  from  3  mm  to  10  mm 
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Table  7:  Comparison  of  LP  performance  for  various  component  thicknesses 


Thickness 

(mm) 

Volume 

(mm^) 

Depth 

(mm) 

Maximum 
Compression  {MPa) 

Maximum 
Tension  {MPa) 

3 

8.8 

0.72 

-555.9 

97.1 

4 

7.3 

0.85 

-576.6 

154.7 

5 

8.1 

0.85 

-596.1 

120.4 

6 

7.2 

0.86 

-604.0 

104.9 

7 

8.8 

0.86 

-587.2 

98.0 

8 

9.1 

0.86 

-620.9 

98.1 

9 

9.8 

0.86 

-623.6 

96.6 

10 

10.1 

0.86 

-623.6 

95.9 

ness  when  the  thickness  is  more  than  4  mm.  The  figure  shows  that  residual 
stress  variations  are  much  higher  for  components  of  3  mm  compared  to  higher 
thickness  components.  The  plots  in  Figure  26  show  the  presence  of  residual 
compressive  stress  at  the  surface  opposite  to  the  peened  surface.  The  magni¬ 
tude  of  these  stresses  is  significant  (>  50  MPa)  when  the  thickness  is  in  the 
range  of  3  to  7  mm.  A  comparison  of  four  performance  measures  is  shown  in 
Table  7.  The  results  indicate  that  initially  (3  to  5  mm)  the  depth  of  compressive 
stress  increases  with  the  thickness  and  then  saturates  upon  further  increases  in 
thickness.  A  potential  explanation  for  these  results  is  that  from  3  to  5  mm, 
the  plastically  affected  depth  is  constrained  by  component  thickness.  From  5 
to  10  mm,  the  plastically  affected  depth  is  limited  by  pressure  pulse  magni¬ 
tude  and  material  properties.  Table  7  shows  that  the  magnitude  of  maximum 
compressive  stress  increases  with  the  component  thickness. 

The  maximum  tensile  stress  tends  to  increase  up  to  4  mm  and  then  starts  to 
decrease  with  increase  in  the  thickness.  Initially,  tensile  stresses  increase  with 
compressive  stresses;  this  is  not  true  after  5  mm  in  thickness.  One  reason 
could  be  that  the  initial  increase  maintains  equilibrium  between  the  tensile 
and  compressive  regions.  Later,  the  volume  of  the  tensile  region  increases  to 
maintain  equilibrium  instead  of  the  magnitude. 

Further  investigations  into  the  results  indicate  that  the  location  of  maximum 
compressive  and  tensile  stress  do  not  change  significantly  with  thickness.  The 
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Figure  27:  Average  residual  stress  field  for  multiple  treatments  of  LP  at  the  same  location 


location  of  maximum  compressive  stress  is  always  on  the  surface  of  the  com¬ 
ponent.  The  location  of  maximum  tensile  stress  is  also  similar  in  all  cases  and 
is  always  below  the  compressive  region.  The  trends  of  compressive  volume  do 
not  seem  to  follow  any  particular  direction.  The  compressive  volume  increases 
when  thickness  varies  from  5  mm  to  10  mm;  however,  it  does  not  follow  any 
specific  trend  in  the  remaining  range  of  thickness.  Although  the  results  ob¬ 
tained  from  the  simulation  indicate  some  trends  and  point  toward  potential 
explanations,  further  investigation  will  be  required. 

3.5.6  Multiple  LP  Treatments  at  the  Same  Location 

A  single  LP  treatment  may  not  be  able  to  achieve  the  depth,  magnitude,  or 
volume  of  compressive  residual  stress  required  for  an  application.  The  depth  of 
compressive  residual  stress  can  be  increased  by  increasing  the  intensity  of  the 
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Figure  28:  Investigated  sequences 


laser  beam,  but  for  certain  applications,  increasing  the  intensity  may  not  be  the 
best  choice  because  of  two  reasons.  After  a  certain  limit,  saturation  may  occur, 
and  the  depth  may  no  longer  be  increased  by  increasing  the  intensity.  Also 
depending  upon  the  material  properties  and  the  process  parameters,  higher 
intensity  can  cause  a  reduction  in  compressive  stress  at  the  center  of  the  spot, 
because  of  the  release  wave.  In  such  cases,  using  multiple  treatments  of  LP  at 
the  same  location  is  an  excellent  option. 

All  simulation  parameters  in  this  study  are  taken  from  convergence  study.  It 
is  believed  that  one  LP  shot  has  the  potential  to  change  the  material  properties 
at  and  around  the  first  laser  spot.  This  effect  can  be  considered  similar  to  a 
welding  process  behavior.  In  welding,  the  material  properties  change  after  one 
weld-pass.  The  physics  of  the  phenomenon  is  approximated  in  the  multi-pass 
welding  simulations.  However,  in  this  research,  for  the  simulation  of  multiple 
treatments  of  LP,  it  is  assumed  that  the  changes  in  material  properties  are  not 
adversely  affecting  the  residual  stress  field  and  therefore  not  considered. 

In  simulations,  the  time  difference  between  two  shots  is  equal  to  the  total  sim- 
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ulation  time.  Here,  the  total  simulation  time  is  approximately  1  to  2  millisec¬ 
onds.  In  real  life  the  time  difference  is  much  higher  and  depends  upon  the 
frequency  of  the  shots.  The  change  in  residual  stresses  is  not  significant  after  1 
millisecond.  The  provided  time  difference  between  two  shots  is  assumed  to  be 
sufficient.  The  simulation  for  the  second  shot  begins  from  the  residual  stress 
distribution  of  the  first  shot.  The  required  properties  of  the  model  are  trans- 
fered  to  the  next  simulation  from  the  previous  simulation.  Figure  27  shows 
the  distribution  of  residual  stress  at  1  mm  from  the  center  line  of  the  spot  after 
three  6.1  GPa  applications  of  LP.  These  results  are  compared  with  published 
experimental  [91  results.  In  the  case  of  the  three  6. 1  GPa  peak  pressure  pulses, 
the  following  results  are  obtained.  The  application  of  the  second  LP  treatment 
at  the  same  location  increases  the  depth  and  the  peak  magnitude  of  the  com¬ 
pressive  residual  stress  by  19.18%  (from  0.73  mm  to  0.87  mm)  and  42.56% 
(from  238.5  MPa  to  340  MPa),  respectively.  The  depth  of  the  compressive 
residual  stress  increases  by  17.24%  (from  0.87  mm  to  1.02  mm)  and  20.9% 
(from  340  MPa  to  409  MPa)  after  the  third  LP  treatment.  In  summary,  im¬ 
provements  obtained  by  using  multiple  treatments  of  LP  at  the  same  location 
are  significant,  but  at  the  same  time,  the  magnitude  of  the  residual  tensile  stress 
also  increases.  Comparison  of  results  for  single  and  multiple  treatments  of  LP 
at  the  same  location  also  shows  that  the  adopted  methodology  is  reasonably 
good  at  predicting  residual  stress  distribution. 

3.5.7  Two  Shots  Sequence  at  Multiple  Locations 

This  section  performs  simulation  of  sequential  LP  at  multiple  locations  for 
two  configurations.  The  first  configuration  considers  two  LP  locations  that 
contact  each  other  without  overlap.  In  the  second  configuration,  there  is  a 
50%  overlap  between  the  two  LP  locations.  The  model  developed  in  section  3 
and  shown  in  Figure  17  is  employed  for  these  investigations.  The  figure  shows 
a  schematic  of  sequential  LP  at  multiple  locations  with  50%  overlap  between 
the  two  locations.  In  the  figure,  the  control  volume  has  two  overlapping  spots 
and  infinite  elements  all  around.  In  the  case  of  0%  overlap,  the  control  volume 
dimensions  are  5  mm,  5  mm,  and  15  mm,  and  in  the  case  of  50%  overlap,  the 
dimensions  are  5  mm,  5  mm,  and  12.5  mm.  Different  sizes  of  control  volumes 
are  used  to  minimize  computational  expenses.  Both  configurations  use  the 
biased  mesh  properties  from  the  mesh  convergence  study. 
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Using  a  Pentium(R)  4  CPU,  3.40  GHz  personal  computer  with  2.00  GB  of 
RAM,  it  takes  approximately  48  Hours  of  CPU  time  to  finish  the  two  shots 
analysis.  Typically,  depending  upon  the  size  and  requirements  of  a  compo¬ 
nent,  it  may  need  anywhere  from  20  to  100  LP  treatments.  Therefore,  the 
time  taken  for  two  shots  is  an  indication  that  the  FE  simulation  for  a  practi¬ 
cal  component  will  be  a  computational  challenge  and  a  suitable  application  of 
surrogate  models  [801  is  required. 

For  sequential  LP  at  multiple  locations,  the  peak  pressure  of  the  pulse  is  5.5 
GPa.  Figure  33  shows  the  residual  stress  distribution  at  the  top  surface  and  at 
a  depth  of  0.5  mm  from  the  surface  of  the  component  for  zero  overlap.  The 
results  shown  in  these  figures  are  averaged.  The  averaging  procedure  is  ex¬ 
plained  in  section  3.  The  figure  shows  that  compressive  stresses  are  reduced  at 
0.5  mm,  compared  to  the  top  surface.  The  As  the  depth  increases,  the  magni¬ 
tude  and  area  of  the  compressive  stress  magnitude  decreases.  This  shows  that 
the  residual  stress  distribution  tends  to  be  uniform  across  the  laser  spot  area, 
it  is  important  to  note  is  that  although  there  is  no  overlap  between  the  two  LP 
locations,  there  seems  to  be  some  interaction  between  the  two  fields.  The  plot 
shows  the  merging  of  the  stress  fields  generated  by  different  shots. 

Figure  34  shows  the  residual  stress  distribution  at  the  surface  and  at  0.5  mm 
from  the  surface  for  50%  overlap.  The  simulation  results  show  that  the  mag¬ 
nitude  of  the  residual  stress  is  greater  in  the  area  of  overlap  as  compared  to 
the  non-overlapping  area.  Compared  to  0%  overlap,  the  magnitude  of  resid¬ 
ual  stress  is  higher.  The  results  indicate  significant  interaction  between  the 
residual  stress  fields  induced  by  both  LP  locations.  It  is  interesting  that  the 
major  effects  of  LP  are  limited  to  a  small  region  and  the  minor  effects  are 
spread  over  the  entire  region.  Investigation  into  the  results  reveal  that  overlap 
between  two  LP  locations  improves  the  magnitude  and  the  depth  of  the  resid¬ 
ual  stress.  Along  with  this  improvement,  overlapping  can  eliminate  the  small, 
shallow  region  of  reduced  compressive  stress  at  the  center  of  the  circular  spot. 
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Figure  29:  Sequence  A  along  the  line 


Figure  30:  Sequence  B  along  the  line 


Figure  31:  Sequence  C  along  the  line 


Figure  32:  Sequence  D  along  the  line 
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Figure  33;  Residual  stress  profile  on  surface  with  0%  Overlap 
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Figure  34;  Residual  stress  profile  on  surface  with  50%  Overlap 
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Figure  35;  Change  in  residual  stress  profile  at  a  line  on  surface  by  seven  shots 
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3.5.8  Seven  Shot  Sequence  at  Multiple  Locations 


In  this  section  seven  shots  are  performed  in  four  different  sequences  to  deter¬ 
mine  the  effect  of  shot  sequence.  A  cylindrical  model  developed  in  section  3 
and  shown  in  Figure  18  is  employed  in  these  investigations.  The  figure  shows 
the  locations  of  seven  shots.  The  four  sequences  investigated  for  comparison 
are  shown  in  different  colors  in  Figure  28.  A  line  is  selected  that  starts  on  the 
surface  and  continuing  into  the  depth,  to  compare  the  residual  stress  profile 
generated  by  each  sequence.  Figures  29  through  32  show  the  residual  stress 
profile  after  each  shot  of  sequences  1  (pink),  2  (light  blue),  3  (blue),  and  4(yel- 
low).  These  plots  show  that  the  effect  of  a  shot  on  the  residual  stress  profile 
depends  upon  the  distance  between  the  point  and  the  shot. 

The  shots  can  be  divided  into  three  categories.  In  the  first  category,  the  resid¬ 
ual  stress  profile  changes  significantly  when  the  point  is  in  the  peened  area. 
These  type  of  shots  tend  to  increase  the  magnitude  of  compressive  stress  as 
well  as  tensile  stress.  In  the  second  category,  changes  in  the  profile  are  mi¬ 
nor  but  positive.  Increases  occur  in  the  magnitude  of  both  compressive  and 
tensile  stresses.  In  the  third  category,  a  shot  tends  to  reduce  the  magnitude  of 
compressive  and  tensile  stresses.  These  reductions  are  very  minor  compared 
to  the  increases  in  the  first  and  the  second  categories.  Further  investigation  is 
required  to  determine  the  exact  distances  of  each  category. 

The  effect  of  multiple  shots  at  the  same  location  are  investigated  in  Section  4.2. 
The  configuration  investigates  only  100%  overlap;  however,  partial  overlap 
(20%,  43%,  75%,  etc.)  is  common  in  many  applications.  The  exploration 
performed  in  this  section  provides  an  oppurtunity  to  investigate  the  effects  of 
number  of  shots  and  partial  overlapping.  Figure  35  plots  the  residual  stress 
profile  along  a  line  on  the  peened  surface  after  each  shot.  The  selected  line 
is  shown  by  a  dotted  line  in  Figure  28.  Figure  35  shows  that  the  first  shot  at 
any  location  induces  the  maximum  magnitude  of  compressive  stress.  In  this 
study,  for  the  combination  of  pressure  pulse  magnitude  (4.5  GPa)  and  material 
(Ti-6A1-4V)  properties,  the  compressive  stress  is  approximately  —450  MPa. 
The  second  shot  at  this  location  increases  the  magnitude  to  —610  MPa.  The 
third  shot  at  the  same  location  with  partial  overlap  increases  it  to  —700  MPa. 

This  results  indicates  that  full  or  partial  overlap  at  a  location  increases  the 
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compressive  stress  magnitude  on  the  surface.  The  increase  after  the  second 
shot  is  higher  than  the  increase  after  third  shot.  Overall  investigation  in  this 
section  shows  that  the  number  of  shots  and  the  layout  of  the  shots  makes  a 
significant  difference  in  the  residual  stress  profile.  However,  sequence  of  shots 
does  not  make  a  considerable  difference  for  the  selected  geometry  and  the 
selected  LP  parameters. 

3.6  Chapter  Summary 

The  section  begins  with  a  standard  formulation  of  a  mixed  variable  optimiza¬ 
tion  problem.  This  is  followed  by  parametric  investigations.  The  parameters 
investigated  are  temporal  variation  of  pressure,  pressure  pulse  magnitude,  spot 
size  and  shape,  component  thickness,  number  of  shots  at  the  same  location, 
overlap,  and  shot  sequence.  The  parameters  were  investigated  with  respect  to 
a  set  of  selected  performance  functions.  The  functions  are  volume  of  compres¬ 
sive  stresses,  maximum  magnitude  of  tensile  and  compressive  stresses,  and 
depth  of  compressive  stresses.  These  parametric  investigations  indicate  that 
our  chosen  performance  functions  have  significant  sensitivities  with  respect  to 
the  LP  parameters.  These  results  motivate  us  to  develop  a  simulation-based 
optimization  strategy  for  obtaining  maximum  benefits  of  laser  peening. 
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4.0  LP  OPTIMIZATION:  ONE  LOCATION 


Optimization  problems  that  result  in  shock,  impact,  and  explosion  type  disci¬ 
plines  typically  have  mixed  design  variables  and  multiple  optimal  solutions.  In 
the  optimization  literature,  many  researchers  have  solved  problems  involving 
mixed  variables  or  multiple  optima,  but  it  is  difficult  to  find  multiple  optimal 
solutions  of  a  mixed- variable  problem.  Therefore,  an  optimization  method  is 
required  that  can  find  multiple  solutions  for  a  mixed- variable  problem. 

A  mixed- variable  niche  particle  swarm  optimization  (MNPSO)  is  developed  to 
solve  such  problems.  The  four  modifications  made  to  the  PSO  are:  Latin  Hy¬ 
percube  sampling-based  particle  generation,  a  mixed- variable  handling  tech¬ 
nique,  a  niching  technique,  and  a  surrogate  model-based  design  space  local¬ 
ization.  The  proposed  method  is  developed  for  a  single  location  laser  peening 
optimization  problem. 

In  many  applications  of  LP,  geometric  configurations  and  dimensional  in¬ 
tegrity  requirements  of  the  component  can  constrain  implementation  of  an  op¬ 
timal  solution.  In  such  cases,  it  is  necessary  to  provide  multiple  alternatives  to 
a  designer  so  that  a  suitable  one  can  be  selected  according  to  the  requirements. 
The  finite  element  analysis  time  of  24  to  72  CPU  hours  for  each  LP  analysis 
makes  it  a  computationally  expensive  optimization  problem  to  solve. 

4.1  Mixed-variable  Niche  Particle  Swarm  Optimization 

Many  engineering  optimization  problems,  such  as  gearboxes,  synchronous 
motors,  and  power  system  operations,  involve  finding  more  than  one  optimal 
solution.  These  optima  can  be  local  or  global.  The  multi-modality  of  the  de¬ 
sign  space  can  cause  difficulties  for  any  optimization  algorithm  being  used  to 
find  optimal  solutions.  This  is  due  to  convergence  toward  local  optima.  In 
many  applications  [81, 821,  it  is  beneficial  to  find  multiple  optimal  solutions. 
The  problem  becomes  difficult  if  the  design  variables  are  a  combination  of 
continuous  and  discrete  variables,  and  the  function  evaluation  is  time  consum¬ 
ing.  The  optimization  method  demanded  by  this  research  should  be  able  solve 
a  multimodal  mixed- variable  optimization  (MMVO)  problem  without  requir¬ 
ing  a  large  number  of  function  evaluations. 
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Based  on  advantages  and  disadvantages  (Section  4.2  section  4)  [83]  of  population- 
based  methods  as  compared  to  gradient-based  methods,  population-based  meth¬ 
ods  are  more  suitable  for  the  problem  at  hand.  Among  population-based  ap¬ 
proaches,  particle  swarm  optimization  (PSO)  [77]  is  selected  over  the  genetic 
algorithm  (GA)  approach  [75,84]  because  of  its  ease  of  implementation,  lower 
number  of  user  parameters,  and  use  of  historical  information. 

The  first  challenge  is  to  incorporate  niching  ability  into  PSO,  and  the  second 
is  to  introduce  a  modification  to  manage  mixed  variables  [85].  Introduction 
of  niching  allows  PSO  to  formulate  and  maintain  subpopulations  converging 
toward  respective  optimum.  Thus,  niching  refers  to  different  optima  within  the 
design  space  that  a  local  subpopulation  explores  and  converges  upon.  There¬ 
fore,  an  optimization  method  with  niching  capabilities  is  able  to  find  and  main¬ 
tain  multiple,  diverse,  and  final  solutions. 

In  the  literature,  many  researchers  have  developed  PSO-based  methods  [83, 
86-90]  to  solve  multimodal  or  mixed- variable  optimization  (MVO)  problems; 
however,  it  is  difficult  to  find  multiple  solutions  of  a  mixed- variable  problem. 
The  GA-based  modified  clearing  method  [91]  demonstrated  niching  using  25 
variables  and  50  peaks  problems.  However,  this  method  can  not  be  used  in  this 
research  because  it  requires  a  large  number  of  function  evaluations.  Parsopolos 
et  al.  [92]  combined  PSO,  neural  network  approximation,  and  threshold  value 
to  incorporate  niching  abilities,  but  the  success  of  this  method  depends  on  the 
threshold  refinement  and  prior  problem  information.  Brits  at  al.  [93]  modified 
this  method,  and  used  particle  tracking  and  subswarm  concept  to  find  multi¬ 
ple  optima.  The  performance  of  this  approach  also  depends  upon  a  threshold 
parameter  and  a  subswarm  radius. 

Seo  et  al.  [87,  90]  developed  the  multigrouped  particle  swarm  optimization 
(MGPSO)  method.  A  traditional  PSO  has  three  terms  to  calculate  particle  ve¬ 
locity,  but  MGPSO  uses  a  fourth  term  to  introduce  niching  properties.  Since 
MGPSO  is  using  an  extra  term  to  introduce  niching  abilities,  its  interference  in 
convergence  is  not  as  much  as  in  other  methods.  Apart  from  the  interference 
with  PSO  convergence,  these  methods  cannot  be  directly  applied  to  mixed- 
variable  optimization  problems. 


84 


To  manage  mixed  variables,  [83]  used  the  two-level  approach.  In  the  two- 
level  approach,  the  mixed  variable  optimization  (MVO)  problem  is  divided 
into  two  problems:  system  level  (mixed- variables)  and  sub- structural  level 
(continuous-variables)  problems.  The  continuous  variable  problem  is  a  sub¬ 
set  of  the  optimization  problem.  Similarly,  [89]  developed  a  hybrid  PSO.  The 
hybrid  PSO  combines  Newton-Raphson  and  PSO;  and  solves  an  MVO  prob¬ 
lem  by  dividing  it  into  a  two-level  problem.  Although  two-level  approach 
shows  that  PSO  can  manage  discrete  variables,  this  approach  cannot  be  im¬ 
plemented  if  the  MVO  problem  cannot  be  easily  divided  into  two  problems. 
Many  researchers  [88,94]  utilized  a  penalty  function  approach  to  handle  mixed 
variables.  This  approach  artificially  modifies  the  design  space,  therefore,  may 
require  a  large  number  of  function  evaluations  for  convergence.  Many  pop¬ 
ulation  based  approaches  simply  round  the  continuous  values  to  their  closest 
discrete  values  to  manage  discrete  variable.  This  approach  can  be  effective  for 
simple  problems  with  wide  range  of  discrete  possibilities.  A  problem  depen¬ 
dent  systematic  approach,  employed  in  this  research,  can  be  more  effective. 

The  literature  survey  shows  that  it  is  possible  to  incorporate  niching  or  mul¬ 
timodal  abilities  into  PSO.  In  this  paper,  PSO  is  modified  to  incorporate  both 
abilities.  An  additional  term  [87, 90]  is  used  in  MGPSO  to  incorporate  nich¬ 
ing  abilities.  This  concept  is  employed  here  because  it  brings  niching  ability 
without  extensively  interfering  in  PSO  convergence.  The  interference  is  lower 
because  the  fourth  term  is  not  always  active.  PSO  convergence  is  similar  to 
typical  PSO  when  additional  term  is  not  active.  In  the  MGPSO  addition  term 
is  not  active  when  a  particle  is  not  intruding  in  a  territory  of  other  group’s 
gbest.  In  the  proposed  method,  the  implementation  of  this  term  is  different. 
The  particles  are  not  divided  into  groups,  and  additional  term  is  active  when  a 
particle  is  in  the  region  of  infiuence  of  a  global  or  local  best  (mbesti)  found  so 
far.  The  mixed  variables  are  managed  by  manipulating  the  random  numbers  in 
the  velocity  calculation  for  the  integer  variable.  Apart  from  these  two  modifi¬ 
cations,  a  third  change  is  introduced:  Latin  Hypercube-  and  surrogate  model- 
based  initial  population  generation.  This  guarantees  that  the  initial  population 
has  particles  from  all  regions  of  the  design  space  and  preserves  the  random 
properties. 

Population-based  methods  require  many  function  evaluations.  Multidisciplinary 
structural  simulation  can  be  computationally  expensive  [95,96].  In  the  litera- 
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ture  [97, 98]  surrogate  models  are  employed  to  mitigate  these  characteristics. 
In  the  proposed  approach,  a  surrogate  model  is  employed  for  an  exhaustive 
search  of  the  design  space.  This  exhaustive  search  indicates  the  location  of 
the  potential  optima.  This  approach  provides  an  effective  starting  population 
without  many  simulations  to  locate  multiple  optimal  solutions.  To  construct 
the  surrogate  model  and  to  perform  multiple  simulations  for  optimization,  a 
parallel  [99, 100]  processing  capability  is  developed.  This  capability  employs 
MATLAB,  PuTTy,  and  the  Ohio  Supercomputing  facilities  to  generate  the  sur¬ 
rogate  model  and  to  run  the  optimization  algorithm.  The  proposed  method  is 
first  validated  on  multiple  test  functions  from  the  literature.  Once  validated, 
the  method  is  implemented  to  solve  the  laser  peening  (LP)  [53, 101, 102]  opti¬ 
mization  problem. 

The  remainder  of  this  section  is  organized  as  follows:  Section  5.2  presents  the 
MMVO  problem  formulation.  The  proposed  MNPSO  is  discussed  in  Section 
5.3.  In  Section  5.4,  the  modified  PSO  is  demonstrated  in  test  problems  and 
results  are  compared  with  the  MGPSO-based  method.  Section  5.5  presents 
a  design  of  experiments  based  surrogate  model  and  the  solution  of  the  LP 
optimization  problem. 

4.2  LP  Multimodal  Problem 

Surface  improvement  by  laser  peening  is  a  trade-off  between  two  opposing  ef¬ 
fects:  the  creation  of  compressive  residual  stress  on  the  surface,  which  tends 
to  increase  life,  countered  in  part  by  the  tensile  stresses,  which  tend  to  de¬ 
crease  life.  Despite  these  two  considerations,  in  many  applications,  such  as 
turbine  blades  and  gears,  preserving  dimensional  integrity  is  essential.  This 
issue  becomes  critical  for  thin  (thickness<  10  mm)  components. 

The  literature  so  far  does  not  apply  to  thin  components  or  to  components  with 
sharp  edges  or  smooth  curved  boundaries.  A  high  LP  pressure  applied  to  a 
thin  target  can  result  in  significant  permanent  deformation  on  the  surface  or 
spall  failure  on  the  back  surface.  A  second  concern  is  that  the  dimensional 
integrity  of  the  component  may  be  at  risk  due  to  high  pressure  magnitude.  The 
use  of  two-sided  peening  to  inhibit  such  losses  is  not  an  option  because  of 
the  interaction  of  the  two  tensile  waves  at  the  center  of  the  component.  This 
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interaction  can  initiate  cracking  and  hinder  the  beneficial  effects.  In  such  cases 
multiple  lower  pressure  shots  are  a  potential  solution.  These  conditions  make 
it  necessary  to  find  multiple  solutions  so  that  multiple  options  can  be  provided 
to  a  designer.  The  designer  can  then  select  a  cost-effective  option  from  the 
solutions. 

4.2.1  LP  Design  Space  for  Multimodal  Optimization 

Residual  stress  field  properties  (volume,  depth,  and  magnitude)  can  be  en¬ 
gineered  by  controlling  LP  parameters  such  as  pressure  pulse  magnitude  and 
shape,  spot  size  and  shape  (continuous),  and  the  number  of  shots  (integer).  The 
five  parameters  that  are  considered  in  optimization  are  pressure  pulse  magni¬ 
tude  (p)  and  midspan  (t),  shot  radius  (r)  and  shape  (.y),  and  number  of  shots 
(n).  Figure  7  shows  the  pressure  pulse  magnitude  and  midspan.  The  pressure 
pulse  magnitude  can  vary  between  2.8  to  8.8  GPa.  The  figure  shows  magni¬ 
tude  on  a  0  to  1  scale.  The  pressure  pulse  shape  is  changed  using  the  pressure 
pulse  midspan.  The  midspan  is  defined  by  the  full  width  at  half  the  maximum 
of  the  pressure  pulse.  The  lower  and  upper  bounds  of  this  variable  (t)  are  30 
and  70  nanoseconds  (ns),  respectively.  For  pressure  pulse  shape,  a  sharp  rise 
to  peak  in  3  ns  is  preselected,  peak  pressure  is  maintained  for  the  next  3  ns, 
and  the  rate  of  drop  is  determined  by  the  variable  (t).  The  range  of  the  variable 
is  such  that  this  can  generate  any  type  of  pulse  between  sudden  (lower  limits) 
and  slow  (upper  limits)  pressure  drops. 

The  size  (r)  and  shape  (^')  of  the  spot  are  other  parameters  that  can  be  con¬ 
trolled.  The  shape  parameter  defined  by  the  diameter  ratio  can  vary  between 
0.5  and  1.0.  The  range  of  radius  is  2.0  to  3.0  mm.  The  combination  of  these 
two  variables  can  generate  many  possible  shapes  and  sizes.  Figure  36  shows 
six  possibilities  of  the  quarter  spot  shape.  The  largest  possible  spot  is  shown 
using  a  black  line  and  a  yellow  triangular  marker.  The  radius  (r)  of  this  shape 
is  3  mm  and  shape  parameter  (s)  is  1.0.  The  smallest  possible  spot  is  shown 
using  black  lines  and  red  circles.  In  this  case,  radius  is  2.0  mm  and  shape  pa¬ 
rameter  is  0.5.  The  number  of  shots  are  indicated  by  n  and  can  vary  between 
1  and  3.  There  are  a  total  of  five  design  variables  {p,r,t,s,n).  To  evaluate 
possible  designs  for  optimization,  finite  element  analysis  of  LP  (Figure  9)  is 
employed.  Due  to  the  small  time  increments  required  to  simulate  the  elastic 
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Figure  36:  Spot  shape  variables 


and  plastic  deformation,  the  computational  cost  of  LP  simulation  is  high. 


4.2.2  Multimodal  Problem  Formulation 

The  goal  of  LP  is  to  generate  as  much  compressive  stress  as  possible  in  the 
component  to  improve  its  surface  properties;  therefore,  the  objective  of  opti¬ 
mization  is  to  maximize  the  volume  of  compressive  stress.  It  is  possible  that 
higher  compressive  magnitude  may  not  be  beneficial  in  all  cases,  but  it  is  al¬ 
most  always  certain  that  the  maximum  magnitude  should  be  on  the  surface. 
To  enforce  this  requirement,  the  compressive  stresses  on  the  surface  are  con¬ 
strained  to  be  at  or  within  a  certain  range  of  the  maximum. 

This  constraint  also  accounts  for  the  nonlinear  behavior  of  LP  for  a  certain 
range  of  parameters,  which  can  generate  a  reduced  compressive  region  at  the 
center  of  the  spot.  This  reduction  may  hinder  the  objective  of  surface  enhance¬ 
ment  techniques.  In  mathematical  terms,  compressive  stresses  at  depth  (d)  are 
constrained  to  be  higher  or  within  a  certain  range  at  the  depth  (d  +  Ad),  where 
(Ad)  is  always  positive.  This  is  a  constraint  on  the  profile  of  the  compressive 
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residual  stresses. 


The  depth  of  compressive  stress  is  also  a  performance  metric  that  should  be  as 
high  as  possible.  The  constraint  that  compressive  stress  depth  should  be  more 
than  1  mm  at  a  critical  location  ensures  this  requirement.  It  is  not  possible  to 
eliminate  tensile  stress,  but  it  can  be  kept  within  a  certain  limit.  The  maximum 
magnitude  of  tensile  stress  induced  by  peening  process  is  controlled  by  the 
third  constraint.  This  constraint  keeps  the  maximum  magnitude  of  the  tensile 
stresses  below  75  MPa. 

The  optimization  problem  statement  is  given  by: 

Maximize  :  Compressive  Volume 

Subjected  to:  Compressive  Stressed)  >  Compressive  Stressed  ^ /^d) 

Compressive  Stress  Depth  >1.0  mm 
Maximum  Tensile  Stress  <75  MPa^ 

Variable  Bounds  :  2.8  <  p  <  8.0  GPa, 

30.0  <t<  70.0  ns, 

0.5  <s<  1.0, 

2.0  <  r  <  3.0  mm, 

1  <  n  <  3, 

This  problem  is  solved  by  using  the  following  optimization  algorithm. 

4.3  Proposed  Method:  Mixed-variable  Niching  PSO 

It  is  difficult  to  find  multiple  optima  by  using  typical  PSO,  because  PSO  has 
an  intrinsic  limitation  that  all  particles  converge  to  only  one  point  at  final  step. 
Even  if  PSO  is  able  to  find  a  local  optimum  in  the  process  of  reaching  a  better 
one,  the  algorithm  is  likely  to  lose  the  local  one.  Therefore,  PSO  needs  to  be 
modified  such  that  it  can  not  only  find  multiple  optima  but  also  retain  the  ones 
it  finds. 

4.3.1  Initial  Particle  Generation 

In  the  first  modification,  the  initial  population  is  generated  using  Latin  Hyper¬ 
cube  Sampling  (LHS)  [1031  instead  of  random  generation.  LHS  is  a  popular 
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method  for  generating  samples.  This  technique  maintains  the  randomness  in 
the  population  and  guarantees  the  coverage  and  exploration  of  the  entire  de¬ 
sign  space.  The  covariance  matrix  required  for  LHS  is  obtained  by  assuming 
that  all  parameters  are  independent  and  normally  distributed.  The  following  is 
a  brief  LHS  procedure  to  generate  a  sample  of  size  p  from  q  variables: 

(i)  Divide  the  design  space  into  p  non-overlapping  intervals  based  on  an 
equal  probability  of  selection  from  each  interval. 

(ii)  Randomly  select  one  value  from  each  interval. 

(hi)  Repeat  Steps  (i)  and  (ii)  for  each  variable  of  q  random  variables  xi ,  X2 , . . . , 
(iv)  Pair  the  p  values  for  each  x/  with  the  p  values  obtained  for  the  other 
at  random. 


4.3.2  Niche  Updating 

The  proposed  method  takes  advantage  of  the  PSO  structure  to  incorporate 
niching  properties,  which  provide  multiple  converged  local  optima.  Two  changes 
are  introduced  to  develop  the  niching  capabilities.  The  first  change  is  a  de¬ 
sign  of  experiments-based  study  to  determine  PSO  parameters  and  the  second 
change  is  an  additional  term  [901  to  traditional  particle  velocity  calculation 
equation. 


Modified  Fourth  Term 

The  traditional  PSO  equation  is: 

=  w/vf  +  earn  {pbesti  -  4)  +  caraigbest  -  4)  (21) 

where  vj  is  the  particle  velocity  at  iteration  k,  Wi  is  the  particle  weight 
function,  c/y  is  the  particle  weight  coefficient  of  each  term  (/^),  r/y  is  a 
random  number  between  0  and  1,  is  the  particle  position  at  iteration  k, 
pbesti  is  the  best  value  of  the  particle,  gbest  is  the  value  of  the  group. 

This  Equation  21  is  modified  to  Equation  22. 
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vf+^  =  Wiv\^Ciirii{pbesti-s\)^Ci2ri2{gbest-s\)^CBri'i{s\-mbesti)  (22) 


where  mbesti  is  the  best  within  the  specified  region  of  influence.  The  new 
constant  c/3  is  calculated  using  Equation  23. 

{1  if  mbesti  —  ^  (l'X\ 

0  Otherwise. 

where  di^mbesu  is  the  distance  between  the  particle  and  mbesti,  d*  is  the  user 
defined  distance  parameter  that  controls  the  region  of  influence. 

To  determine  mbesti,  all  the  particles  are  sorted  in  descending  order  of  their 
objective  function  values.  In  the  sorted  particles  the  first  particle  is  the  mbesti. 
Any  particles  that  is  within  d*  distance  of  mbesti  can  not  be  mbesti.  In  the 
sorted  particles,  the  next  highest  particles  that  is  not  in  the  region  of  influence 
of  mbesti  is  the  mbesti-  Similarly,  other  mbesti  and  their  region  of  influence 
are  determined. 

The  fourth  term  concept  is  taken  from  MGPSO  paper  but  the  implementation 
of  the  term  is  different.  In  MGPSO  the  term  is  active  when  a  particles  intrudes 
the  territory  of  the  other  group’s  gbest.  In  the  proposed  method  the  term  is 
active  only  when  more  than  q  particles  are  within  the  distance  {d*)  of  the 
mbesti.  The  distance  {d*)  and  the  number  of  particles  q  are  parameters  the 
user  must  provide  at  the  beginning  of  the  algorithm.  These  parameters  can 
be  easily  calculated  using  the  problem  information.  For  example,  for  a  two- 
dimensional,  50-peak  function,  the  number  of  particles  is  200,  d*  is  0.05,  and 
p  is  4. 

If  the  particle  is  in  a  region  of  influence,  its  velocity  calculation  includes  the 
fourth  term;  otherwise  it  does  not.  A  particles  is  considered  to  be  in  the  region 
of  influence  of  mbesq,  if  the  distance  between  the  particle  and  mbesq  is  less 
than  d*.  The  particles  that  are  moved  based  on  this  criterion  are  expected  to 
improve  exploration  of  the  unexplored  regions  of  the  search  space.  Based  on 
the  generated  random  number  r/3,  the  particle  is  relocated  outside  the  region. 
This  strategy  allows  PSO  particle  to  converge  on  the  nearest  optimum  and  pre- 
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serve  multiple  optimal  solutions.  This  relocation  may  try  to  send  the  particles 
out  of  the  design  space.  To  prevent  this  situation,  particles  are  relocated  inside 
the  design  space. 


Cognitive  and  Social  Coefficients 

The  two  constants  c/i  and  c/2  play  a  significant  role  in  preserving  the  optima 
already  found.  A  design  of  experiments-based  parametric  study  is  performed 
to  determine  the  combination  of  these  two  parameters  that  is  best  suited  to 
preserve  these  optima.  In  the  study,  the  constants  c/i  and  c/2  is  changed  from 
0.0  to  2.5  in  the  step  of  0. 1.  This  is  a  full  factorial  design.  There  are  two  factors 
and  26  levels.  The  constant  values  are  further  refined.  This  study  provides  an 
effective  parameter  setting  for  a  problem.  This  study  provides  a  range  for 
both  constants  instead  of  a  value.  Any  value  of  the  constants  within  this  range 
can  preserve  the  multiple  optimal  solutions.  For  a  two-dimensional  10-peak 
function,  approximately  1.2  (c/i)  and  0.01  (c/2)  are  the  best  combination  for 
preserving  the  multiple  optimal  solutions. 

For  \D  and  2D  problems  the  constants  c/i  and  c/2  have  the  same  value  through 
all  iterations.  For  higher  dimensional,  the  constant  values  are  changed  with  the 
iteration  number  [104,  1.  The  constants  are  changed  linearly.  The  cognitive 

f 

factor  (c/i)  starts  with  a  low  value  and  linearly  increases  to  The  social 

f 

factor  (c/2)  starts  at  a  high  value  c^2  linearly  decreases  to  a  low  value  c-2- 
The  parameters  are  updated  according  to 

(24) 

^max 

Here,  kmax  is  the  maximum  number  of  iterations  and  k  is  the  current  iteration. 

The  upper  and  lower  bound  of  the  two  constants  are  adjusted  for  each  function 

f  f 

independently.  For  a  3D,  50  peaks  problem,  the  constants  ck ,  ck ,  c%,  and 
are  1.5,  2.0,  0.1,  and  0.01. 

The  parameter  setting  provided  by  the  study  gives  no  or  minimal  importance 
to  the  social  behavior  of  the  particles  (0.01  for  ID  and  between  0.1  and  0.01 
for  3D).  Typically,  the  strength  of  the  population-based  approach  comes  from 
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Iteration  History 


Figure  37;  Number  of  Peaks  Found  w.r.t.  Iteration  Number 


its  social  behavior.  The  purpose  of  the  social  behavior  is  to  explore  the  design 
space.  In  case  of  multimodal  optimization  problems,  it  is  necessary  to  explore 
the  design  space  but  it  is  also  important  to  converge  to  local  optima.  The  low 
value  of  the  social  constant  reduces  the  exploration.  This  reduction  does  not 
affect  the  performance  of  the  algorithm  significantly  because  of  the  additional 
term  in  the  PSO  equation  and  the  parameters  q.  The  additional  term  and  the 
parameters  q  forces  particles  to  unexplored  design  space.  This  argument  can 
be  verified  by  monitoring  the  number  of  peaks  found  at  a  iteration. 

A  3Z)  problem  with  50  peaks  is  solved.  The  problem  details  and  the  PSO 
parameters  are  shown  in  Section  4.2.  Figure  37  shows  the  number  of  peak 
found  at  each  iteration.  The  figure  shows  that  number  of  peaks  found  are 
increasing  with  the  number  of  iterations.  One  of  the  reason  for  finding  more 
peaks  as  iteration  increases  is  that  modified  PSO  is  forcing  particles  to  explore 
unexplored  design  space.  This  result  shows  that  the  modifications  are  bringing 
niche  ability  with  minor  effects  on  the  exploration  capabilities  of  the  PSO. 

Advantages  and  Disadvantages 

This  method  has  multiple  advantages  and  disadvantages  over  the  multigrouped 
particle  swarm  optimization  (MGPSO)  method.  The  first  advantage  over  MG- 
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PSO  [87]  is  that  in  the  proposed  method  particles  are  not  divided  into  groups. 
Therefore  the  group  parameters  and  calculations  required  for  to  form  are  not 
needed.  The  second  advantage  is  that  there  is  a  limit  to  the  maximum  number 
of  particles  {q)  that  can  converge  to  an  optimal  solution.  This  forces  more  par¬ 
ticles  to  explore  the  design  space.  In  MGPSO,  there  is  no  limit  on  maximum 
number  of  particles  in  a  group.  A  disadvantage  is  that  the  proposed  method 
needs  an  additional  parameter  (maximum  number  of  particles  around  a  peak, 
q).  The  calculation  of  this  parameter  is  very  simple.  This  is  approximated 
by  dividing  the  total  number  of  particles  by  expected  number  of  optima.  The 
third  advantage  is  that  the  proposed  method  uses  dynamic  coefficients  instead 
of  constant.  This  approach  modifies  the  particle  behavior  depending  upon  the 
iteration  number.  This  modification  of  constants  assist  the  particles  to  con¬ 
verge  to  multiple  optima. 

4.3.3  Integer  Variable  Technique 

A  minor  alteration  in  the  random  numbers  r/y  of  the  velocity  calculation  (Equa¬ 
tion  22)  can  help  manage  the  integer  variables.  Due  to  the  spacing  and  range 
of  the  integer  variable,  there  are  only  a  limited  number  of  possibilities.  Based 
on  the  random  numbers,  one  possibility  is  selected.  For  example,  if  there  are 
two  possibilities  for  a  random  number,  below  or  equal  0.5  indicates  to  the  first, 
and  above  0.5  indicates  to  the  second  possibility.  If  there  are  three  possibilities 
the  threshold  will  be  0.33,  if  four  possibilities  then  0.25,  instead  of  0.5.  This 
modification  allows  both  continuous  and  integer  variables  to  be  used  in  the 
procedure  with  no  inconsistency. 
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k+l 


+  Ciirn  {pbesti  -  4)  +  camigbest  -  4)  +  c/3r/3 (4  -  mbesti) 
0.5  X  (1)  +  1  X  ra{2  -  1)  +  1  X  r/2(3  -  1)  +  1  x  ro(l  -  2) 

0.5  +  1  X  {0.0, 0.5, 1.0}  X  (2  -  1)  +  1  X  r/2(3  -  1)  +  1  x  r/3(l  -  2) 
0.5  +  1  X  0.5  X  (1)  +  1  X  rail)  +  1  x  rB(-l) 

0.5  +  0.5  +  1  X  {0.0,0.25,0.50,0.75,1.0}  x  (2)  +  1  x  r/3(-l) 

0.5  +  0.5  +  1  x0.25x  (2)  +  l  xr/3(-l) 

0.5  +  0.5  +  0.5  +  0.5x(-l) 

0.5  +  0.5  +  0.5-0.5 
1 


The  integer  handling  technique  is  explained  using  a  example.  It  is  assumed 
that  at  iteration  k,  Wi  =  0.5,  =  I,  =  1,  pbesti  =  2,  gbest  =  3,  mbesti  =  2, 

Oil  =  1,  c/2  =  1  and  c/3  =  1.  These  values  are  substituted  in  equation  22  and 
the  velocity  at  iteration  ^  +  1  is  calculated.  As  shown  in  the  example,  three 
random  numbers  (r/i,  ra,  and  r/3)  are  needed  to  calculate  the  velocity  . 
The  random  numbers  are  decided  in  the  sequence  of  r/i,  ra,  and  then  r/3.  There 
are  three  options  0.0,  0.5,  and  1.0  for  r/i.  These  options  will  make  the  second 
term  0.0,  0.5,  and  1.0,  respectively.  Let  us  assume  that  a  random  number 
generation  selected  0.5  out  of  three  possibilities.  Therefore,  the  second  term 
is  0.5.  Similarly,  there  are  five  options  for  ra.  For  ra,  out  of  five  possibilities 
random  number  generation  selected  0.25.  This  makes  the  third  term  0.5.  In 
the  fourth  term,  the  random  number  r/3  has  only  one  option.  The  option  is  0.5. 
The  velocity  is  1.  Similar  procedure  is  used  to  handle  integer  variables. 

4.3.4  Parallel  Processing 

PSO  and  surrogate  model  development  requires  multiple  finite  element  sim¬ 
ulations  of  the  system  under  consideration.  For  example,  one  simulation  can 
take  anywhere  from  24  to  72  hours  of  CPU  time  using  the  Glenn  cluster  at 
OSC.  To  manage  this  cost,  an  algorithm  is  developed  (using  MATLAB  and 
PuTTy)  that  can  submit  multiple  simulations  in  parallel  to  OSC  and  the  extract 
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Figure  38:  Multimodal  Test  Function  1  Figure  39:  Multimodal  Test  Function  2 


results.  This  insures  that  computational  and  licensing  resources  are  optimally 
utilized.  A  detailed  description  of  the  process  is  given  in  section  4  of  the  report. 


4.4  Multimodal  Test  Problems 

Before  solving  the  problems,  it  is  important  to  keep  in  mind  that  all  the  pro¬ 
posed  modifications  are  not  required  for  each  function  considered  in  this  re¬ 
search.  Depending  upon  the  function,  one,  two,  or  all  three  modifications  can 
be  incorporated. 

4.4.1  A  Periodic  Function  with  Peaks  of  Equal  Size  and  Interval 

The  first  multimodal  test  function  is  a  \D  problem  with  five  peaks  of  equal 
height  at  equal  interval.  The  mathematical  formula  for  the  function  is  shown 
in  Equation  25  and  the  function  is  plotted  in  Figure  38.  The  variable  range  is 
0  to  1  and  peaks  are  at  0.1,  0.3,  0.5,  0.7,  and  0.9. 

P\{x)  =  sin^  {5tix)  (25) 

The  modification  to  manage  integer  variable  is  not  used  here  because  the 
variable  is  a  continuous  one.  The  constants  d*,  g,  number  of  particles,  and 
numbers  of  iterations  used  in  solving  this  problem  are  0.08,  3,  20,  and  100, 
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respectively.  All  the  parameters  are  selected  after  investigating  effects  of  each 
parameter.  For  example,  the  number  of  iterations  is  selected  to  be  100  because 
this  makes  sure  that  the  algorithm  is  able  to  maintain  the  discovered  optima  for 
a  large  number  of  iterations.  The  number  of  particles  is  20  which  is  sufficient 
to  solve  the  problem.  The  circles  in  Figure  38  show  the  final  population  of  the 
modified  PSO.  The  figure  shows  that  the  method  is  not  only  able  to  find,  but 
also  to  maintain  the  optimal  solutions. 

4.4.2  A  Periodic  Function  With  Peaks  of  Unequal  Size  and  Interval 


Figure  40;  Multimodal  Bump  Function  Figure  41;  Particles  for  50  Peak  Function 

This  multimodal  test  function  has  five  peaks  of  unequal  height  at  unequal 
intervals.  The  function  is  shown  in  Figure  39.  The  mathematical  formula  of  the 
function  is  given  in  Equation  26.  The  locations  of  the  five  peaks  are  0.0796, 
0.246,  0.4491,  0.6786,  and  0.929  and  the  corresponding  function  values  are 
0.9895,  0.886,  0.6576,  0.3782,  and  0.1585,  respectively. 

P2(x)  =  g-2(ln2)(^)2  (5^[_^0.75  _  ^^26) 

The  constants  d*,  g,  number  of  particles,  and  number  of  iterations  used  in  solv¬ 
ing  this  problem  are  0.08,  3,  20,  and  100,  respectively.  The  circles  in  Figure  39 
show  the  particles  after  100  iterations  of  modified  PSO  on  the  function.  The 
proposed  method  can  find  and  maintain  all  the  optima  for  a  specified  number 
of  iterations.  This  problem  is  solved  again  after  changing  the  number  of  itera¬ 
tions  to  1000  and  keeping  all  the  other  parameters  same.  The  final  distribution 
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of  the  particles  after  1000  iterations  is  similar  to  the  distribution  after  100  it¬ 
erations.  Similar  tests  have  been  performed  on  all  numerical  problems.  These 
tests  prove  that  the  number  of  iterations  does  not  matter  for  the  algorithm  in 
maintaining  all  the  optimum  solutions. 

4.4.3  Multimodal  Bump  Function 

The  previous  two  test  functions  are  one-dimensional  functions  with  only  five 
peaks.  Bump  function  is  a  parametric  multimodal  function,  in  which  the 
dimension,  number  of  peaks,  radius  of  peaks,  curvature  of  peaks,  height  of 
peaks,  and  location  of  peaks  are  parametric.  This  variability  allows  us  to  gen¬ 
erate  many  multimodal  functions.  The  mathematical  formula  for  the  bump 
function  is  given  as 


where  is  the  Euclidean  distance  of  the  particle  from  the  midpoint  of 
peak,  rk,  h^,  and  are  the  radius,  height  and  shape  parameter,  respectively  of 
the  peak,  and  K  is  the  number  of  peaks.  With  the  bump  function,  multiple 
problems  with  different  levels  of  difficulty  can  be  generated.  The  difficulty 
of  the  multimodal  problem  can  be  changed  using  parameters  K,  h^,  r^,  and 
dk-  The  top  view  of  a  two-dimensional  bump  function  is  shown  in  Figure  40. 
The  number  of  peaks  K  is  20.  The  range  of  each  variable  is  [0, 1].  The  center 
of  each  peak  is  randomly  initialized,  with  only  one  constraint,  the  distance 
between  the  midpoints  of  any  two  and  peaks  is  greater  than  or  equal  to 
in  +  rm).  The  figure  shows  a  function  with  random  location  Xk,  height  h^,  and 
shape  parameter 

The  results  of  the  presented  method  and  the  MGPSO  method  are  compared 
[871.  For  comparison,  2,  3,  and  5  dimensional  bump  functions  with  10,  20,  30, 
40,  and  50  peaks  are  generated.  The  height  of  the  peak  is  a  random  number 
between  0.8  to  1.0  for  all  the  functions.  The  peak  curvatures  range  between 
1.5  to  2.0  for  all  the  functions.  The  radius  is  different  for  each  function.  De¬ 
pending  upon  the  number  of  peaks  and  dimension  of  the  function,  a  radius 
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Table  8:  Radius  Range  for  Bump  Functions 


Radius  ranges  of  peaks 

Dim. 

10 

20 

30 

40 

50 

2D 

0.08  to  0.10 

0.07  to  0.09 

0.06  to  0.08 

0.055  to  0.075 

0.05  to  0.07 

3D 

0.27  to  0.30 

0.19  to  0.21 

0.16  to  0.18 

0.14  to  0.15 

0.13  to  0.15 

5D 

0.48  to  0.53 

0.40  to  0.44 

0.32  to  0.35 

0.30  to  0.33 

0.28  to  0.31 

Table  9:  Number  of  Particles  to  Solve  the  Function 


Number  of  Particles 

Dim. 

10 

20 

30 

40 

50 

2D 

50 

75 

100 

120 

150 

3D 

60 

80 

250 

500 

600 

5D 

600 

1000 

1200 

1500 

2000 

range  is  selected.  The  selected  ranges  are  shown  in  Table  8.  The  radius  of  a 
peak  is  a  random  number  in  the  range. 

These  generated  functions  are  solved  using  the  proposed  method  and  MGPSO- 
based  algorithm.  The  multimodal  functions  and  the  PSO  parameters  are  the 
same  for  both  methods  except  the  parameters  that  are  not  common  in  both 
methods.  In  the  case  of  the  proposed  method,  the  third  term  is  applied  to  a 
selected  group  of  particles,  whereas  in  the  case  of  the  competing  method  it  is 
applied  to  all  the  particles.  The  selection  criteria  are  defined  with  the  method 
description. 

Parametric  studies  have  been  performed  for  each  function  to  find  the  regions  of 
influence  (d*)  and  p,  and  the  number  of  particles.  Table  9  shows  the  number 
of  particles  used  in  each  of  the  functions.  The  number  of  particles  for  each 
function  is  selected  using  a  manual  parametric  investigation.  Figure  41  shows 
the  final  population  of  the  modified  PSO  on  the  2D  test  function  with  50  peaks. 
The  figure  demonstrates  that  modified  PSO  particles  are  converging  to  all  the 
optima,  are  able  to  find  optima,  and  are  able  to  maintain  it  for  100  iterations. 
To  demonstrate  the  robustness  of  the  algorithm,  each  problem  is  solved  10 
times  with  a  different  initial  population  of  PSO. 


99 


Table  10:  Bump  Function:  Method  with  MGPSO  third  term 


Number  of  Peaks  Found 

10 

20 

30 

40 

50 

Dim. 

SR 

SD 

SR 

SD 

SR 

SD 

SR 

SD 

SR 

SD 

2D 

8.8 

0.75 

18.8 

0.75 

28.0 

1.41 

33.2 

1.17 

41.2 

1.78 

3D 

7.6 

1.11 

13.0 

1.18 

25.9 

1.45 

34.5 

1.96 

40.4 

1.63 

5D 

7.0 

1.10 

12.2 

1.89 

14.2 

2.44 

20.2 

3.25 

24.8 

3.28 

The  results  the  of  competing  and  proposed  methods  are  shown  in  Tables  10 
and  11,  respectively.  The  shown  results  are  the  average  success  rate  of  the  10 
mns.  The  average  success  rate  and  the  standard  deviation  of  the  success  rate 
is  reported.  In  the  tables,  Dim.,  SR,  and  SD  indicate  the  function  dimension, 
success  rate,  and  standard  deviation  in  the  success  rate,  respectively.  The  suc¬ 
cess  rate  is  defined  by  the  number  of  optima  found  by  the  optimizer.  The  SD 
is  shown  to  demonstrate  that  the  optimizer  is  consistent  in  its  performance.  An 
optimum  is  considered  to  be  found  if  at  least  one  particle  is  within  a  certain 
distance  of  the  peak.  This  distance  is  selected  to  be  15%  of  the  peak  radius. 
These  results  demonstrate  that  the  method  is  able  to  find  multiple  optima  con¬ 
sistently. 

A  comparison  of  the  results  shows  that  both  methods  are  able  to  solve  the 
multimodal  optimization  problems.  The  performance  of  both  methods  is  sim¬ 
ilar  on  two-dimensional  problems.  However,  for  three-dimensional  problems, 
the  proposed  method  performs  slightly  better  than  the  MGPSO-based  method. 
The  results  of  the  five-dimensional  problem  shows  that  the  proposed  method 
finds  a  higher  number  of  optima  compared  to  the  other  method.  The  compar¬ 
ison  of  standard  deviation  in  the  success  rate  shows  that  the  proposed  method 
is  more  consistent  in  finding  multiple  solutions. 


4.5  Multimodal  LP  Optimization  and  Results 

The  LP  parameters,  mesh  properties,  material  model,  and  geometric  boundary 
conditions  of  the  3D  model  can  be  found  in  section  3.  In  brief,  the  LP  simu¬ 
lation  employed  in  multimodal  optimization  is  a  symmetric  3D  model  shown 
in  Figure  9.  This  simulation  considers  pressure  pulse  magnitude,  duration. 
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Table  1 1 :  Bump  Function:  Modified  PSO  with  boolean  third  term 


Number  of  Peaks  Pound 

10 

20 

30 

40 

50 

Dim. 

SR 

SD 

SR 

SD 

SR 

SD 

SR 

SD 

SR 

SD 

2D 

9.6 

0.66 

19.6 

0.66 

29.9 

0.30 

39.6 

0.66 

49.9 

0.30 

3D 

9.5 

0.67 

19.2 

0.60 

28.6 

0.92 

36.5 

1.20 

43.6 

1.43 

5D 

9.1 

0.70 

18.5 

0.92 

25.3 

0.78 

33.2 

1.25 

39.4 

1.11 

and  shape,  spot  shape  and  size,  and  number  shots  at  the  same  location  as  op¬ 
timization  variables.  To  manage  the  simulation  time,  a  design  of  experiments 
(DOE)-based  surrogate  model  is  employed. 


4.5.1  Design  of  Experiments-based  Surrogate  Model 

In  many  traditional  methods,  approximation  is  constructed  around  a  current 
design  point  to  obtain  the  next  point.  However,  in  the  case  of  LP,  only  a  limited 
number  of  simulations  can  be  performed;  therefore,  a  surrogate  model  that 
covers  the  all  of  the  move  limits  is  required  for  each  step.  Traditionally,  a 
model  is  required  to  predict  objective  functions  and  constraints.  However,  in 
the  case  of  LP,  the  residual  stress  values  at  all  the  nodes  need  to  compute  the 
volume  and  profile  of  the  compressive  stress,  magnitude  of  the  tensile  stress, 
and  depth  of  the  compressive  stress.  Therefore,  a  surrogate  model  is  needed 
for  each  node  response.  In  the  3D  model,  there  are  68,921  (41  nodes  per  5 
mm)  or  531,441  (81  nodes  per  5  mm)  nodes,  and  constructing  the  response 
surface  for  residual  stress  at  each  node  is  necessary  to  calculate  performance 
metrics  (objective  and  constraints).  This  surrogate  model  combines  DOE  and 
a  multi-point  response  surface. 

To  develop  a  surrogate  model,  a  3-level  full  factorial  design  is  considered. 
Simulations  are  performed  to  calculate  residual  stress  values  at  all  the  nodes 
in  all  designs  of  DOE.  Since  there  are  five  variables,  a  total  of  (3^  =)  243 
simulations  are  performed.  This  set  up  provides  the  matrix  [X]  (Eq.  29)  and 
the  response  vector  [Y]  for  each  node  to  construct  a  response  surface.  Eor  each 
node  a  surrogate  model  is  assumed  as  shown  in  the  following  equation: 
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^  —  Po  +  Pl/l  W  +  •••  +  ^kfk{^)  +£ 


(28) 


where  P/,  /  =  0, 1, 2,  are  the  unknown  coefficients,  8  is  residue  in  the  re¬ 
gression  model,  and  fi{xj)  are  the  functions  of  individual  or  combined  x/.  The 
least  squares  method  is  used  to  find  unknown  coefficients,  and  it  yields 


(29) 

f  =  xp,  e  =  Y-Y  (30) 

The  above  computation  (Eqs.  29  and  30)  is  performed  for  each  node  response. 
These  calculations  provide  unknown  coefficients  P  and  Y  for  all  nodes.The  first 
advantage  of  this  surrogate  model  is  that  this  provides  the  residual  stress  field 
of  the  entire  model  instead  of  an  objective  or  constraint  value.  This  field  can 
be  used  to  calculate  any  quantity  needed  and  can  also  be  used  for  the  purpose 
of  fatigue  life  calculation.  The  second  advantage  is  that  once  the  model  is 
developed,  it  does  not  need  updates  at  every  point  unless  algorithm  reaches  a 
converged  point  and  wish  to  change  or  reduce  move  limits.  This  is  because  the 
surrogate  model  covers  the  entire  range  of  design  variables. 

4.5.2  Results  and  Discussion 

All  three  modifications  of  PSO  are  employed  to  solve  the  LP  optimization 
problem.  The  value  of  d*  is  1.0,  and  the  value  of  p  is  2.  The  value  of  p  is  2 
because  the  number  of  PSO  particles  is  20.  The  number  of  particles  is  limited 
to  20  based  on  available  resources.  The  variable  ranges  of  pressure  magnitude 
and  the  number  of  shots  drive  the  selection  of  d*.  To  incorporate  compres¬ 
sive  stress  depth  and  profile  and  tensile  stress  constraints,  a  penalty  function 
approach  is  employed.  The  penalty  is  proportionate  to  the  magnitude  of  the 
constraint  violation.  The  penalty  is  such  that  a  feasible  solution  has  a  better  ob¬ 
jective  function  compared  to  an  infeasible  solution.  This  is  a  modified  version 
of  exterior  penalty  function  used  in  population-based  optimization  methods. 

Two  iterations  of  the  modified  PSO  are  used.  The  first  iteration  uses  the  surro¬ 
gate  model.  This  iteration  is  a  type  of  exhaustive  search  of  the  design  space. 
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Table  12:  LP  Optimization  Results 


Sol.  # 

P 

n 

r 

t 

s 

Obj. 

Tmax 

Cmax 

Depth 

S 

1 

3.62 

3 

2.95 

38.1 

0.63 

4.78 

74.7 

-612.0 

1.50 

100.9 

2 

5.35 

2 

2.22 

46.8 

0.69 

4.93 

74.5 

-652.7 

1.35 

132.5 

3 

5.59 

2 

2.78 

51.8 

0.52 

4.97 

74.0 

-703.7 

1.25 

148.3 

4 

6.00 

1 

2.76 

65.5 

0.93 

5.56 

74.9 

-603.7 

1.25 

119.5 

because  a  function  evaluation  does  not  require  significant  computational  time. 
This  iteration  identifies  the  potential  solutions/regions  in  the  design  space, 
which  become  part  of  the  initial  population  for  the  second  iteration.  For  the 
second  iteration,  these  potential  solutions  (particles)  are  combined  with  the 
particles  generated  using  an  LHS  technique.  The  second  step  begins  with  a 
population  that  has  representatives  of  potential  optimum  solutions.  Full-scale 
simulations  are  used  in  the  second  iteration.  Table  12  shows  the  optimization 
results  provided  by  the  modified  PSO  algorithm. 

In  the  table,  Obj.  indicates  compressive  volume,  Tmax  indicates  maximum  uni¬ 
directional  tensile  stress,  Cmax  indicates  maximum  uni-direction  compressive 
stress,  and  S  indicates  the  maximum  principal  stress.  The  algorithm  is  able  to 
find  multiple  solutions.  The  first,  second,  third,  and  fourth  solutions  require  3, 
2,  2,  and  1  peening  shots,  respectively.  The  first  solution  uses  3.62  GPa,  the 
second  uses  5.35  GPa,  the  third  uses  5.58  GPa,  and  the  fourth  uses  6.0  GPa 
pressure  pulse  magnitudes. 

The  pressure  pulse  magnitude  and  the  number  of  shots  are  the  major  param¬ 
eters  that  vary  in  the  solutions.  There  is  a  direct  relationship  between  the 
pressure  magnitude  and  the  number  of  shots:  higher  pressure  magnitudes  re¬ 
quire  lower  number  of  shots  to  produce  similar  effects.  All  the  solutions  satisfy 
the  profile,  tensile  stress,  and  compressive  stress  depth  constraints.  The  active 
constraints  are  the  maximum  tensile  stress  and  the  profile  of  the  compressive 
stress.  The  maximum  tensile  stress  constraint  is  active  for  all  the  solutions. 
The  profile  constraint  is  active  for  the  second,  third,  and  fourth  solutions.  A 
closer  examination  of  the  results  provides  the  following  information: 

(i)  There  are  multiple  solutions  for  the  formulated  LP  problem. 

(ii)  Multiple  shots  at  lower  pressure  at  the  same  location  create  a  higher  depth 
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of  compressive  stresses  compared  to  single  shot  of  higher  pressure. 

(hi)  One  shot  at  higher  pressure  creates  higher  magnitude  compressive  stress 
compared  to  multiple  shots  at  lower  pressure. 

(iv)  An  elliptical  spot  shape,  when  compared  to  a  circular  spot  shape,  induces 
a  lower  magnitude  of  tensile  stress  for  the  same  compressive  stress  mag¬ 
nitude. 

A  non-niching  optimization  method  will  find  only  the  last  solution  because  it 
has  the  highest  objective  function  value.  The  other  solutions  obtained  using 
modified  PSO  can  potentially  be  more  useful  for  certain  applications.  The  ob¬ 
jective  function  value  of  these  solutions  is  not  the  same.  This  shows  that  first 
three  solutions  are  local  optima.  However  with  these  solutions,  a  designer  has 
multiple  options  to  choose  from.  Depending  upon  the  effects  of  the  LP  param¬ 
eters  on  a  component,  an  appropriate  solution  can  be  selected.  For  example, 
when  peening  leading  edge  of  a  turbine  blade,  where  minor  dimensional  vari¬ 
ations  can  affect  its  performance  significantly,  the  first  solution  may  have  an 
advantage  over  the  last  solution.  For  an  application  where  the  LP  parameters’ 
effects  are  within  the  allowable  tolerances,  the  last  solution  can  be  selected. 
If  a  high  magnitude  of  compressive  stress  is  required,  then  the  third  solution 
is  better  than  other  solutions.  If  both  depth  and  magnitude  of  compressive 
stress  are  equally  important,  then  the  second  solution  is  better  than  the  other 
solutions. 

Another  advantage  of  multiple  solutions  is  that  different  equipment  can  be 
used  to  accomplish  a  set  goal.  For  example  if  a  company  has  multiple  ma¬ 
chines  with  different  capacities  (laser  intensity /pressure  magnitude).  One  of 
the  machine  can  be  busy  due  its  capabilities  and  not  possible  to  use  for  all 
jobs.  In  such  cases,  a  solution  can  be  picked  based  on  the  available  machine 
capabilities.  This  allows  the  company  to  optimally  utilize  the  available  re¬ 
sources  due  to  the  multiple  solutions.  Therefore,  each  of  these  solutions  have 
advantages  and  disadvantages  over  the  other  solutions,  and  depending  upon 
the  requirements  the  designer  can  select  the  one  that  is  appropriate. 
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4.6  Chapter  Summary  Remarks 


A  MNPSO  is  proposed  to  solve  the  multimodal  mixed-variable  optimization 
problem.  The  MNPSO  employs  LHS-based  population  generation,  a  niching 
strategy,  and  a  mixed- variable  handling  technique.  LHS-based  generation  is 
used  to  obtain  a  diverse  population,  the  niching  strategy  is  used  to  find  and 
preserve  multiple  optima,  the  mixed-variable  technique  and  surrogate  model 
design  space  localization  is  used  to  manage  continuous  and  integer  variables. 
An  exhaustive  search  is  performed  using  the  surrogate  model  to  localize  poten¬ 
tial  multiple  solutions.  The  method  is  demonstrated  on  multiple  test  problems 
and  on  laser  peening  optimization.  The  solution  of  the  laser  peening  problem 
provides  a  designer  with  multiple  options.  The  designer  can  select  a  suitable 
option  depending  upon  dimensional  integrity  requirements.  The  successful 
demonstration  on  multiple  problems  shows  that  the  proposed  method  can  be 
employed  to  solve  multimodal  mixed- variable  problems. 
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5.0  LP  OPTIMIZATION:  MULTIPLE  LOCATIONS 


In  this  research,  a  progressive  simulation-based  design  optimization  strategy 
is  developed  that  can  be  applied  to  highly  nonlinear  impulse-type  processes 
such  as  shot  peening,  laser  peening,  and  bullet  impacts  on  aircraft  structural 
components.  The  design  problems  entail  the  use  of  multiple  fidelities  in  sim¬ 
ulation,  time-consuming  elastic-plastic  analysis,  and  mixed  types  of  optimiza¬ 
tion  variables.  An  optimization  strategy  based  on  progressively  increasing  the 
complexity  and  fidelity  is  developed,  along  with  suitable  surrogate  models. 
Multi-level  fidelity  models  include  axi-symmetric,  symmetric  3D,  and  full- 
scale  simulations  to  enable  design  optimization.  The  first  two  models  are  used 
to  perform  parametric  studies  and  to  localize  the  potential  design  space.  This 
creates  a  reduced  design  space  and  an  effective  starting  point  for  the  subse¬ 
quent  optimization  iterations.  All  steps  employ  the  modified  particle  swarm 
optimization  for  mixed  variables. 

The  design  methodology  is  demonstrated  on  LP  of  a  structural  component. 
The  LP  parameters,  pressure  pulse  and  spot  dimensions,  impulse  locations  (all 
continuous),  number  of  shots  (integer),  and  the  location  of  shots  (discrete)-are 
the  optimization  variables  with  stress  constraints. 

5.1  Progressive  Multifidelity  Optimization  Strategy  (PMOS) 

Multidisciplinary  optimization  of  structures  subjected  to  high-speed  impact 
processes  is  complex  due  to  the  nonlinear  and  transient  nature  of  the  finite  el¬ 
ement  analysis.  The  simulations  consist  of  elastic-plastic  analysis  with  small 
time  steps  and  large-scale  finite  element  models  for  practical  aircraft  struc¬ 
tures.  The  experiments-based  design  approaches  traditionally  used  are  cost 
prohibitive  for  parameter  optimization.  Therefore,  a  simulation-based  design 
methodology  is  required.  In  this  research,  a  design  methodology  is  developed 
for  structures  subjected  to  the  shock-type  loadings  that  occur  in  shot  peen¬ 
ing,  laser  peening,  explosions,  and  bullet  impacts.  The  optimization  problem 
is  also  challenging  due  to  the  presence  of  the  continuous  and  discrete  design 
variables  [851,  a  mixed- variable  optimization  (MVO)  problem. 

In  the  literature,  depending  upon  the  geometry,  loads,  and  boundary  condi- 
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tions,  multi-fidelity  simulation  and  surrogate  models  are  used  to  solve  aircraft 
structural  analysis  problems.  Livne  et  al.  [1051  used  polynomial-based  equiva¬ 
lent  plate  modeling  techniques  for  wing  structures  analysis.  This  work  demon¬ 
strates  that  2D  models  can  be  used  to  approximate  the  nonlinear  aeroelastic 
response  of  a  wing-box  subjected  to  in-plane  compressive  forces.  Robinson  et 
al.  [1061  used  physics-based  and  finite  element  models  with  parameter  map¬ 
ping  on  a  wing  design  problem.  This  research  proved  that  low-fidelity  and 
high-fidelity  models  connected  with  parameter  mapping  can  assist  an  opti¬ 
mization  approach  and  save  computational  time.  Chen  et  al.  [1071  employed 
multiple  objective-oriented  surrogate  models,  each  constructed  for  a  local  re¬ 
gion  and  combined  using  a  boolean  operations.  This  approach  updates  the 
surrogate  models  based  on  the  preexisting  surrogate  model  and  confidence  in¬ 
terval  on  approximation. 

Glaz  et  al.  [1081  compared  a  weighted  average  of  multiple  surrogate  mod¬ 
els  (design  of  experiments,  polynomial  response,  kriging,  radial  basis  neural 
network)  to  each  individual  surrogate  model  and  demonstrated  the  advantage 
of  the  weighted-average  surrogate  [109, 1101.  Although  these  techniques  are 
effective,  they  may  not  be  the  best  options  for  impulse-type  mixed-variable 
problems.  In  this  research,  a  progressive  multifidelity  optimization  strategy 
(PMOS)  is  developed.  This  strategy  combines,  low-  and  high-fidelity  simula¬ 
tions  [105, 1061  and  respective  surrogate  models  [107, 1081,  as  well  as  mixed- 
variable  strategy  [851  to  solve  impulse-type  problems. 

The  basis  of  PMOS  is  that,  depending  upon  the  problem  information  and  as¬ 
sumptions,  multiple  simulation  models  can  be  developed  to  approximate  the 
impulse  response.  Multiple  impulses  may  be  required  in  a  process,  but  the 
effects  of  a  single  impulse  at  a  location  can  be  determined  using  a  single  axi- 
symmetric  model  to  address  a  large-scale  problem.  However,  assumptions  in 
the  model  limit  the  number  of  variables  that  can  be  considered  in  the  model. 
Similarly,  in  a  symmetric  2D  model,  the  number  of  assumptions  can  be  re¬ 
duced,  so  this  model  can  incorporate  additional  variables.  Simulation  of  an  ar¬ 
bitrary  structural  component  with  many  impulses  at  multiple  locations  requires 
an  analysis  model  without  symmetric  assumptions.  Typically,  such  models  are 
computationally  expensive  but  can  consider  most  of  the  variables  in  the  prob¬ 
lem.  These  qualities  of  the  models  are  utilized  in  developing  the  optimization 
strategy.  The  uniqueness  of  the  proposed  method  is  that  it  localizes  the  design 
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Figure  42:  Multifidelity  optimization  strategy 

space  using  axi- symmetric  2D  and  symmetric  3D  simulation  and  surrogate 
models,  and  determines  the  optimum  solution. 

As  shown  in  Figure  42,  the  PMOS  involves  a  three-step  hierarchical  procedure. 
The  design  space  is  localized  by  parametric  studies  and  by  solving  optimiza¬ 
tion  problems  with  low  fidelity  models,  which  also  assist  in  eliminating  insen¬ 
sitive  optimization  parameters.  These  steps  provide  an  effective  starting  point 
and  significantly  reduce  the  search  space  for  the  third  step  of  the  optimization 
procedure.  The  reduced  design  space  requires  fewer  function  evaluations  in 
the  subsequent  optimization  step.  When  employed  together  these  three  steps 
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reduce  the  number  of  the  full  3Z)  model  simulations. 


The  developed  method  is  for  impulse-type  problems.  There  will  always  be  an 
integer  variable  in  such  problems,  for  representing  the  number  of  impulses.  As 
shown  in  section  5,  the  PSO  [77,831  is  modified  to  generate  diverse  initial  pop¬ 
ulation  and  handle  mixed  variables.  The  modifications  introduced  are  Latin 
Hypercube  Sampling-based  initial  particle  generation  and  a  discrete  variable 
handling  technique. 

In  the  first  modification,  the  initial  population  is  generated  Latin  Hypercube 
Sampling  (LHS)  instead  of  random  generation.  This  modification  is  explained 
in  Section  5.3.1  of  section  5.  The  second  modification  is  to  handle  integer  vari¬ 
able.  This  modification  is  described  in  section  5.3.3  of  section  5.  Since  there 
is  not  an  inbuilt  termination  criterion  for  PSO,  depending  upon  the  objective 
requirements  and  available  resources,  the  maximum  number  of  function  eval¬ 
uations  terminates  the  optimization. 

5.1.1  Sub-parametric  Surrogate  Model 

An  important  property  of  impulse-type  problems  is  the  repetition  of  impulses. 
A  surrogate  model  is  developed  for  the  3D  model  that  utilizes  this  repetitive¬ 
ness.  This  surrogate  model  is  based  on  the  characteristic  of  an  impulse-type 
problem  that  each  impulse  has  certain  set  of  parameters  that  define  it.  If  two 
impulse  locations  on  a  structure  are  significantly  apart  from  each  other  and  the 
parameters  defining  these  impulses  are  the  same,  then  both  impulses  generate 
similar  local  effects.  The  proposed  model  exploits  this  quality  of  the  impulse 
response.  The  major  advantage  of  this  idea  is  that  it  does  not  require  simula¬ 
tion  of  all  impulses  at  different  locations  to  approximate  the  cumulative  effect 
of  all  impulses. 

For  a  traditional  design  of  experiments-based  surrogate  model,  a  certain  num¬ 
ber  of  simulations  are  performed.  The  results  of  these  simulations  are  used 
to  construct  an  approximation.  This  is  not  the  best  approach  for  nonlinear 
impulse-type  problems  due  to  the  computational  cost  of  each  simulation.  In 
this  surrogate  model,  instead  of  performing  a  complete  simulation  for  each 
impulse,  simulations  for  two  or  three  impulses  are  performed  to  extract  the  ef¬ 
fects  of  individual  parameters.  A  three-impulse  simulation  requires  less  time 
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Configuration  1 


Figure  43:  Two  possible  layouts 


than  a  five-impulse  simulation  while  still  providing  the  information  required 
to  build  a  database  for  the  surrogate  model.  Once  all  the  individual  parameters 
and  interaction  effects  are  extracted,  they  can  be  combined  to  approximate  the 
performance  for  most  configurations. 

The  above  concept  is  explained  using  two  examples.  The  first  configuration  in 
Figure  43  shows  a  three-impulse  layout.  The  effects  for  this  layout  can  be  ap¬ 
proximated  using  simulation  of  only  two  impulses,  because  the  third  impulse 
does  not  directly  overlap  the  first  impulse.  Similarly,  the  second  configura¬ 
tion  in  Figure  43  has  a  five-impulse  layout.  Since  the  maximum  number  of 
directly  interacting  impulses  is  three,  multiple  three-impulse  simulations  can 
be  used  to  approximate  the  resulting  performance.  The  surrogate  models  re¬ 
quires  shots  1,  2,  and  4  to  be  simulation  using  FEA.  The  residual  stress  profile 
for  shot  3  can  be  approximated  from  FE  simulation  of  shot  2.  Similarly,  the 
residual  stress  profile  for  shot  5  can  be  approximated  from  FE  simulation  of 
shot  4.  After  the  database  is  generated,  it  can  be  used  to  approximate  solution 
for  the  required  layout. 

Although  this  surrogate  model  reduces  the  computational  cost,  it  has  two  dis¬ 
advantages.  The  first  disadvantage  is  that  it  partially  ignores  the  effects  of  im¬ 
pulse  sequence.  The  second  disadvantage  is  that,  irrespective  of  the  distance 
between  the  two  impulses  unless  they  overlap,  the  interaction  effect  is  ne- 
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glected.  These  drawbacks  are  not  significant  when  compared  with  the  advan¬ 
tages  of  the  surrogate  model.  Moreover,  since  full  simulations  are  performed 
after  each  optimization  iteration  convergence,  the  surrogate  model  serves  as  an 
intermediary  interactive  tool  and  the  final  results  are  always  reached  without 
making  any  of  the  above  approximations. 

5.1.2  PMOS:  Advantages  and  Disadvantages 

The  proposed  strategies  have  certain  advantages  and  disadvantages.  The  first 
advantage  is  that  the  proposed  strategy  utilizes  axi-symmetric  2D  and  sym¬ 
metric  3D  models.  The  second  advantage  is  that  the  computational  cost  is 
lower  than  the  full  3D  models.  The  third  advantage  is  that  the  stepwise  pro¬ 
cess  identifies  infeasible  regions  without  performing  full  3D  simulation.  A 
disadvantage  is  that  it  is  possible  to  lose  some  interaction  effects  among  some 
parameters.  This  disadvantage  can  be  avoided  by  performing  parametric  in¬ 
vestigations  before  attempting  the  design  optimization.  Based  on  these  in¬ 
vestigations,  the  optimization  strategy  can  be  modified  to  accommodate  the 
interaction  effects. 

The  fourth  advantage  is  that  mixed  variables  are  managed  effectively  without 
inconsistency.  The  fifth  advantage  is  that  the  proposed  surrogate  model  does 
not  require  simulation  of  all  the  impulses.  Overall,  the  strategy  solves  an  MVO 
problem  at  a  lower  computational  than  a  traditional  approach,  without  sacri¬ 
ficing  effectiveness.  The  strategy  is  demonstrated  in  the  laser  peening  (LP) 
process. 

5.2  Multiple  Locations  Peening  Optimization 

The  pressure  pulse  magnitude  and  shape,  spot  size  and  shape,  spot  location, 
and  number  of  shots  are  the  optimization  variables.  There  are  multiple  chal¬ 
lenges  in  implementing  the  three-step  strategy  to  optimize  laser  peening  with 
these  variable.  The  first  challenge  is  to  develop  a  simulation  model  that  takes 
the  required  LP  parameters  into  consideration.  A  parametric  model  is  required 
that  can  not  only  consider  the  above  variables  but  also  the  location  of  the 
shots,  the  amount  of  overlap,  and  the  sequence  of  the  shots  (mixed  variables). 
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The  following  sections  present  models  available  in  the  literature  and  develop 
a  model  to  implement  the  optimization  strategy. 

5.2.1  Multiple  Simulation  Models 

Finite  element  simulations  of  a  structure  experiencing  impulses  can  be  per¬ 
formed  using  multiple  models  with  respective  assumptions.  Depending  upon 
the  geometry  and  boundary  conditions,  response  to  an  impulse  can  be  mod¬ 
eled  using  axi- symmetric,  symmetric  3Z),  or  full-scale  3D  models.  Typically, 
axi-symmetric  and  symmetric  3Z)  models  require  lower  computational  times 
than  full  3Z)  models.  The  full-scale  3D  model  requires  a  higher  computational 
time  but  can  consider  more  optimization  variables  than  the  other  two  models. 
A  detailed  description  of  axi-symmetric  2D  and  symmetric  3D  models  can  be 
found  in  section  3. 

3D  Model 

The  2D  and  symmetric  3D  simulations  have  applications  in  one  or  multiple 
shots  at  the  same  location.  Because  a  typical  application  may  require  10  to 
100  locations  or  shots  on  a  large  surface  (50  to  500  mm),  an  FE  model  of  such 
a  system  can  be  computationally  prohibitive  [1021.  To  investigate  the  effects 
of  shots  at  more  than  one  location,  the  different  overlapping  configurations, 
and  different  sequences  of  shots,  further  improvements  are  needed  in  the  FE 
model.  A  plate  model  is  developed  to  investigate  required  parameters.  The 
model  is  shown  in  Figure  19.  The  details  of  the  model  are  in  Section  3.7  of 
section  3. 

5.3  Problem  Formulation,  PMOS  Implementation,  Validation,  and  Re¬ 
sults 

5.3.1  Optimization  Formulation:  MVO 

For  impulse-type  loading,  the  objective  of  an  optimization  problem  can  be  to 
minimize  cost  of  the  process,  minimize  or  maximize  damage,  or  maximize 
a  performance  function.  The  objective  for  LP  [6,  36,  101,  111,  1121.  is  to 
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maximize  the  compressive  stress  volume  in  the  peened  component.  The  first 
constraint  forces  the  compressive  stresses  on  the  surface  to  be  at  or  within 
a  certain  range  of  the  maximum.  This  constraint  is  included  to  mitigate  the 
effects  of  the  nonlinear  behavior  of  the  process,  which  within  certain  range 
of  parameters  can  generate  a  reduced  compressive  region  at  the  center  of  the 
spot.  This  reduction  may  hinder  the  objective  of  the  surface  enhancement 
techniques. 

In  mathematical  terms,  compressive  stresses  at  depth  (d)  are  constrained  to  be 
higher  or  within  a  certain  range  at  the  depth  (d  +  Ad),  where  (Ad)  is  always 
positive.  The  second  constraint  requires  that  the  depth  of  compressive  stress 
is  more  than  1  mm  at  a  selected  location.  The  maximum  magnitude  of  tensile 
stress  induced  by  the  peening  process  is  controlled  by  the  third  constraint. 
This  constraint  keeps  the  maximum  magnitude  of  the  tensile  stresses  below 
150  MPa.  The  design  variables  are  pressure  pulse  magnitude  (p/)  and  shape 
{t\i  and  t\i),  spot  shape  (si)  and  size  (r/),  location  of  shot  (x/  and  y/),  and  number 
of  shots  at  the  same  location(n/).  Here  i  indicates  the  shot.  The  lower  and 
upper  bounds  for  each  variable  are  shown  in  the  formulation: 


Maximize  :  Compressive  Stress  Volume 

Subjected  to:  Compressive  Stress{(P  >  Compressive  Stress{d Ad) 

Compressive  Stress  Depth  >1.0  mm 
Maximum  Tensile  Stress  <150  MPa, 

Variable  Bounds  :  2.8  <  p/  <  8.0  GPa, 

15.0  <  tii  <  50.0  ns, 

50.0  <  t2i  <  150.0  ns, 

0.5  <  Si  <  1.0, 

2.0  <  r/  <  3.0  mm, 

2.5  <  Xi  <  7.5  mm, 

2.5  <  Vi  <  7.5  mm, 

I  <ni<  3 

A  traditional  optimization  can  solve  the  problem  formulated  above.  This  will 
use  full  scale  five  shot  plate  simulations.  A  plate  simulation  requires  4  CPU 
days  on  the  Glenn  cluster  at  OSC.  A  simulation  for  three  shots  at  each  of  the 
five  locations  requires  12  days  of  CPU  time.  Due  to  the  computational  time 
issues,  finding  an  optimal  solution  using  the  proposed  optimization  strategy  is 
more  effective  than  traditional  optimization. 
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5.3.2  Optimization  Strategies  Implementation 


The  three  models  of  LP  required  for  the  PMOS,  axi- symmetric,  symmetric 
3D,  and  the  parametric  plate-require  approximately  20,  170,  and  5760  CPU 
minutes,  respectively.  As  a  starting  point  for  the  three-step  optimization  strat¬ 
egy,  parametric  investigations  [1131  are  performed  to  determine  the  effects  of 
individual  parameters.  These  investigations  find  a  significant  variation  in  para¬ 
metric  sensitivities  with  respect  to  performance  metrics.  In  the  first  step  of  the 
strategy,  four  LP  parameters  are  optimized  using  2D  simulation.  The  second 
step  considers  four  LP  parameters  (pressure  pulse  magnitude  and  duration, 
spot  shape,  and  number  of  shots)  and  uses  symmetric  3D  simulation.  In  the 
third  step,  four  parameters  (pressure  pulse  magnitude,  spot  location  and  size, 
number  of  shots  at  the  same  location)  are  considered. 

5.3.3  Step  1:  Optimization  using  a  2D  Model 

In  this  section,  the  first  optimization  step  is  implemented  employing  2D  simu¬ 
lation.  There  are  three  LP  parameters  in  the  optimization:  pressure  magnitude 
(p),  shot  radius  (r),  and  pressure  pulse  shape  (t).  As  shown  in  Figure  44,  the 
pressure  pulse  shape  is  defined  using  two  design  variables  (ti,t2);  as  a  result, 
there  are  a  total  of  four  design  variables  (p,  r,t\,t2).  Based  on  the  capabilities 
of  the  equipment,  the  peak  pressure  (p)  range  is  3.5  to  8.8  GPa,  and  the  spot 
radius  range  is  2.0  to  3.0  mm.  For  pressure  pulse  shape,  a  sharp  rise  (peak  in 
3  ns)  is  preselected,  and  the  rate  of  drop  is  determined  by  the  variables 
The  ranges  of  the  variables  are  such  that  these  can  generate  any  type  of  pulse 
between  a  sudden  (lower  limits)  and  slow  (upper  limits)  pressure  drop.  The 
ranges  of  the  two  shape  variables  are  15  to  50  ns  {t\)  and  90  to  150  ns  (u). 
The  optimization  problem  statement  is  given  by: 
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Figure  44:  Optimization  Formulation 


Maximize  :  Compressive  Stress  Volume 

Subjected  to:  Compressive  Stress{c^  >  Compressive  Stress{d Ad) 

Compressive  Stress  Depth  >  1.0  mm 
Maximum  Tensile  Stress  <  150  MPa, 

Variable  Bounds  :  3.5  <  /?  <  8.8  GPa, 

2.0  <  r  <  3.0  mm, 

15.0  <ti<  50.0  ns, 

90.0  <t2<  150.0  ns. 

The  number  of  particles  and  generations  in  PSO  are  set  at  20.  The  values  c/i 
and  c/2  are  taken  as  1.7.  The  optimization  results  show  that  a  peak  pressure 
of  5.42  GPa,  a  spot  radius  of  2.1  mm,  and  pressure  shape  parameters  of  17.0 
ns  and  150.0  ns  are  needed  to  achieve  the  desired  residual  stress  profile.  The 
formulation  can  achieve  the  required  depth  and  limit  the  maximum  tensile 
stresses  within  the  specified  limit.  The  maximum  tensile  stress  constraint  is 
the  active  constraint  at  the  optimum.  It  can  be  seen  from  the  results  that  the 
low  spot  radius  can  keep  the  reduced  compressive  region  at  the  center  of  the 
spot  to  within  the  limit.  The  first  shape  parameter  (ti)  is  near  the  lower  bound 
of  the  variable  range,  and  the  second  shape  parameter  (^2)  is  at  the  upper  bound 
of  the  variable  range.  This  indicates  a  sharp  drop,  but  a  long  pressure  pulse  is 
favored  for  the  formulation. 
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5.3.4  Step  2:  Optimization  using  a  symmetric  3D  model 

In  this  section  the  second  step  of  the  strategy  is  implemented.  Pressure  pulse 
magnitude  (p/),  number  of  shots  (n),  and  laser  shot  shape  (si)  are  the  variables. 
The  pressure  pulse  shape  is  taken  from  the  previous  step,  and  the  pressure 
magnitude  range  is  reduced  for  better  exploration  of  the  design  space.  This 
model  obtains  a  effective  starting  point  from  the  previous  step  but  experiences 
computational  challenges.  A  surrogate  model  is  used  to  reduce  the  computa¬ 
tional  expense.  A  detailed  description  of  the  model  can  be  found  in  section 
5.5.2  in  section  5. 


Optimization  Formulation  and  Results 

In  pressure  pulse  magnitude  (p/)  and  spot  shape  (si)  the  subscript  i  denotes 
the  shot  of  the  process.  In  the  following  formulation,  the  volume  of  the 
compressive  region  is  the  objective  function.  All  the  constraints  and  variable 
bounds  for  fatigue  strength  are  shown  in  the  equation: 

Maximize  :  Compressive  Stress  Volume 

Subjected  to:  Compressive  Stress{d)  >  Compressive  Stress{d -\- Ad) 

Compressive  Stress  Depth  >1.0  mm 

Maximum  Tensile  Stress  <150  MPa, 

Variable  Bounds  :  4.0  <  p/  <  6.5  GPa, 

0.5  <  Si  <  1.0, 

2.0  <  r  <  2.5  mm, 

1  <  <  3, 

The  results  of  the  problem  show  that  two  shots  with  pressure  pulse  magnitudes 
of  3.6  and  4.7  GPa,  radius  of  2.28  mm,  and  shape  parameters  of  1.0  and  0.80 
are  the  optima.  The  pressure  magnitude,  radius,  and  shape  variables  are  guided 
by  the  profile  constraint.  In  this  optimization  the  profile  constraint  is  an  active 
constraint. 
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Table  13:  Results  for  methodology  validation 


Problem 

Model 

P 

r 

ti 

Number 

used 

{GPa) 

(mm) 

{ns) 

{ns) 

Problem  1 

2D 

5.72 

2.49 

33.5 

150.0 

Problem  1 

3D 

5.58 

2.44 

38.1 

110.1 

Problem  2 

2D 

8.00 

2.29 

33.1 

150.0 

Problem  2 

3D 

8.00 

2.39 

26.6 

149.1 

5.3.5  Methodology  Validation 

An  assumption  of  the  proposed  optimization  strategy  is  that  a  comparatively 
lower  fidelity  model  {2D)  can  provide  a  solution  similar  to  that  of  higher  fi¬ 
delity  models  {2D  and  plate).  To  validate  this  assumption,  the  problem  for¬ 
mulated  for  the  2D  model  is  solved  using  both  2D  and  2D  models.  The  same 
PSO  parameters  are  used  for  both  models.  The  termination  criterion  in  both 
approaches  is  the  maximum  number  of  iterations.  This  is  because  a  popula¬ 
tion  based  technique  can  easily  determine  the  potential  region  but  can  require 
a  large  number  of  iterations  to  find  the  converged  optima. 

The  results  from  both  models  are  shown  in  Table  13.  The  solutions  in  the 
table  show  that  results  from  both  models  are  similar.  The  table  shows  that 
the  parameters  pressure  magnitude  (p=  5.72  and  5.58  GPa),  and  radius  (r  = 
2.49  and  2.44  mm)  are  similar  for  both  models.  The  pressure  pulse  shape 
parameters  {ti  and  ti)  are  different  in  both  solutions.  In  the  2D  solution  ti 
is  greater,  and  in  the  2D  solution  is  greater.  These  differences  cancel  each 
other  out  to  some  extent,  creating  similar  results. 

A  second  problem  is  also  solved  using  both  models.  The  second  problem  is 
the  same  as  that  of  the  first  problem,  except  that  the  constraint  on  the  residual 
stress  profile  is  removed  (the  first  constraint  from  the  formulation  in  Section 
3.4).  Table  13  shows  that  the  results  from  both  models  are  similar.  The  re¬ 
sults  of  the  second  problem  indicate  that  the  pressure  magnitude  has  a  direct 
relationship  with  the  compressive  stress  profile.  The  pressure  magnitude  value 
is  at  the  maximum  possible.  The  radius  and  shape  are  not  at  maximum  to 
control  the  reduced  compression  at  the  center  of  the  spot.  This  study  shows 
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that  2D  simulation  can  be  employed  to  save  computational  time  and  that  the 
assumption  made  in  the  optimization  strategy  is  reasonable. 

5.3.6  Step  3:  Optimization  using  a  Parametric  Plate  Model 

At  this  point,  the  first  two  steps  of  the  optimization  strategy  have  been  imple¬ 
mented.  The  third  step  of  the  strategy  is  important  because  the  model  in  this 
step  considers  the  greatest  number  of  variables  of  any  model.  For  example  in 
LP,  it  is  necessary  to  include  the  impulse  location  variables  that  can  generate 
various  overlapping  configurations  or  layouts.  To  do  that,  it  is  necessary  to  use 
the  parametric  plate  model.  However,  the  required  4-day  CPU  time  forces  us 
to  use  the  proposed  surrogate  model. 


Surrogate  Model 

The  parameters  that  define  an  LP  shot  are  pressure  pulse  magnitude,  shape,  and 
duration,  spot  shape  and  size,  and  amount  of  overlap.  Two  LP  shots  with  the 
same  parameters  at  two  locations  generate  similar  local  residual  stress  fields. 
Therefore,  the  proposed  surrogate  model  can  be  applied  to  the  LP  problem. 
For  the  parametric  plate,  three  shot  simulations  are  performed  to  extract  the 
effects  of  individual  parameters.  Information  collection  take  less  time  using 
three  shot  simulations  than  five- shot  simulations.  Once  the  individual  parame¬ 
ters  and  interaction  effects  are  extracted,  they  can  be  combined  to  approximate 
the  residual  stress  profile  for  most  configurations. 

An  example  of  the  surrogate  model  is  shown  in  Figure  45.  This  example  is  a 
schematic  for  only  one  variable  (pressure  magnitude).  Suppose  a  design  needs 
a  residual  stress  field  for  6.0  GPa  pressure  magnitude.  The  database  is  avail¬ 
able  for  4.5,  5.5,  6.5,  7.5,  and  8.5  GPa  pressure  magnitudes.  The  database 
for  each  magnitude  contains  residual  stress  values  at  the  plate  nodes.  To  ap¬ 
proximate  the  residual  stress  field  for  6.0  GPa  ,  the  residual  stress  fields  from 
5.5  and  6.5  are  used.  Equation  31  shows  the  formulation  to  approximate  the 
residual  stress  field.  Here  7?^^)  indicate  the  residual  stress  field  at  pressure  p. 
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Figure  45:  Schematic  of  combining  database  information  to  obtain  desired  stress  field 


Similarly  an  example  for  two  variables,  the  five  simulations  at  4.5,  5.5,  6.5, 
7.5,  and  8.5  GFa  are  performed,  and  the  local  residual  stress  field  at  the  laser 
spot  is  analyzed  and  saved.  For  the  amount  of  overlap,  simulations  are  per¬ 
formed  for  0,  25,  50,  75  and  100%  overlapping  configurations.  The  effects 
of  the  second  and  the  third  shots  at  the  same  location  are  different  from  each 
other;  therefore,  the  100%  overlapping  configration  for  the  third  shot  is  per¬ 
formed  seprately.  To  generate  the  residual  stress  field  for  five  shots  on  different 
locations  with  20%  overlap  for  2”^  to  5^^  shots  and  6. 1  GPa  pressure  for  all 
shots,  the  following  process  is  used.  The  effects  of  the  first  shot  are  directly  in¬ 
terpolated  from  the  local  residual  stress  fields  of  5.5  and  6.5  GFa  magnitudes. 
The  second  shot  is  interpolated  from  the  local  residual  stress  fields  of  5.5  and 
6.5  GFa  pressure  pulses  and  the  local  residual  stress  fields  of  the  0  and  25% 
overlap  configurations. 

Similarly,  the  approximate  residual  stress  field  can  be  generated  for  any  num¬ 
ber  of  shots  on  the  plate.  In  this  research,  linear  interpolation  is  used;  however, 
depending  upon  the  parametric  investigations,  this  can  be  changed  to  higher- 
order  interpolations.  This  approach  does  not  need  updates  within  the  move 
limits.  Just  as  with  the  surrogate  model  for  the  full  3D  model,  this  also  pro¬ 
vides  the  residual  stress  field  for  all  the  nodes  instead  of  a  selected  node.  A 
limited  number  of  five- shot  simulations  are  needed  to  construct  the  database 
because  the  fifth  shot  is  often  located  at  the  center  of  the  plate  overlapping  the 
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Table  14:  Comparison  of  results  from  Surrogate  and  FEA  models 


1 

2 

3 

4 

FEA 

App. 

FEA 

App. 

FEA 

App. 

FEA 

App. 

cv 

34.7 

26.7 

43.1 

32.8 

44.9 

32.4 

42.5 

33.5 

MC 

-419.7 

-479.8 

-563.1 

-612.1 

-529.6 

-623.6 

-479.2 

-398.3 

CD 

0.63 

0.54 

0.63 

0.63 

1.00 

1.13 

0.75 

0.75 

MT 

77.1 

93.8 

97.6 

98.5 

84.8 

85.3 

110.0 

97.0 

previous  four  shots. 


Surrogate  Model  Validation 

The  following  investigations  were  performed  to  compare  the  estimates  of  the 
surrogate  model  with  FE  simulations.  Four  random  peening  configurations 
are  selected.  FE  analysis  is  performed  to  determine  the  compressive  stress 
volume,  depth,  and  maximum  magnitude.  The  surrogate  model  is  also  used  to 
determine  these  quantities.  Here,  pi  indicates  the  pressure  magnitude  of  the 
shot,  Xi  indicates  the  x  co-ordinate  of  the  shot  location,  and  y/  indicates 
the  y  co-ordinate  of  the  shot  location.  The  peening  parameters  in  the  first 
configuration  are  p\  =  4.1,  p2  =  5.3,  p^  =  4.6,  p^  =  5.1,  p^  =  3.9,  xi  =  3.6, 
X2  =  5.8,  X3  =  2.9,  X4  =  5.6,  X5  =  4.1,  yi  =  3.8,  y2  =  3.2,  y3  =  6.4,  y4  =  5.8, 
and  y5  =  6.0.  The  peening  parameters  in  the  second  configuration  are  4.7,  5.6, 
4.3,  5.6,  4.7,  3.3,  6.4,  4.0,  5.6,  5.2,  3.1,  3.4,  5.7,  6.5,  and  5.3.  The  peening 
parameters  in  the  third  configuration  are  5.6,  4.2,  4.8,  4.6,  6.0,  2.6,  6.5,  4.3, 
5.7,  5.0,  3.4,  2.9,  5.8,  5.9,  and  4.6.  The  peening  parameters  in  the  fourth 
configuration  are  5.5,  4.1,  5.7,  5.1,  5.1,  2.7,  6.4,  2.5,  6.4,  4.1,  4.4,  3.1,  5.9, 
6.2,  and  4.0. 

Table  14  shows  the  results  from  both  approaches.  The  comparison  shows  that 
the  approximation  and  the  FE  models  do  not  match  exactly.  However,  the 
surrogate  model  is  able  to  capture  the  trends.  The  minimum  errors  in  the 
compressive  stress  volume,  depth,  and  maximum  magnitude  and  tensile  stress 
maximum  magnitude  are  21.2%,  0.0%,  8.7%,  and  0.06%.  The  maximum  er¬ 
rors  are  27.8%,  14.3%,  17.7%,  and  21.1%,  respectively.  Overall,  the  proposed 
surrogate  model  has  the  potential  to  be  used  in  the  optimization.  An  important 
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Figure  46;  Optimal  layout  of  the  peening  process 


aspect  to  note  is  that  the  finite  element  estimation  took  4  CPU  days,  while  the 
surrogate  model  estimation  took  just  10  CPU  minutes. 


Optimization  Formulation  and  Results 

The  optimization  formulation  and  results  of  the  third  step  are  provided  below. 
Similar  to  the  results  of  the  previous  two  steps,  the  results  show  that  a  lower 
pressure  magnitude  is  favored  because  of  the  residual  stress  profile  constraint. 
The  spot  layout  is  shown  in  Figure  46.  In  this  problem,  the  constraint  on  the 
residual  stress  profile  and  the  maximum  tensile  stresses  restricts  the  maximum 
compressive  stress  magnitude  and  the  compressive  stress  volume: 

Maximize  :  Compressive  Stress  Volume 

Subjected  to:  Compressive  Stressed)  >  Compressive  Stress{d h,d) 
Compressive  Stress  Depth  >1.0  mm 
Maximum  Tensile  Stress  <  150  MPa, 
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Table  15;  Results  of  optimization  strategy 


Shot  No. 

Pressure 
magnitude  (GPa) 

Location 
X  (mm) 

Location 
y  (mm) 

Radius 

(mm) 

No.  of 

shots 

1 

4.42 

2.5 

2.87 

2.28 

1 

2 

4.67 

6.5 

2.5 

2.28 

1 

3 

4.53 

2.72 

6.38 

2.28 

1 

4 

6.20 

6.28 

6.5 

2.28 

1 

5 

4.37 

5.76 

5.08 

2.28 

1 

Variable  Bounds  :  2.5  <  x/  <  7.5  mm, 

2.5  <  3^/  <  7.5  mm, 

3.5  <  Pi  <  6.5  GPa, 

1  <  <  3, 

The  optimization  results  are  given  in  Table  15.  In  addition  to  achieving  the  set 
objective  and  complying  with  the  constraints,  the  optimization  is  able  to  find 
a  layout  that  covers  the  region  of  interest.  The  region  of  interest  is  in  between 
2.0  and  7.5  on  both  (x  and  y)  axes  and  is  shown  by  a  red  dotted  line  in  Figure 
46.  The  results  indicate  that  for  the  formulated  problem,  a  pressure  magnitude 
of  approximately  4.5  GPa  can  achieve  the  objective.  The  magnitude  of  the  4^^ 
shot  is  different  from  the  other  shots  at  6.2  GPa.  A  possible  reason  for  this  is 
that  the  5^^  shot  overlaps  the  4^^  shot.  Because  of  this  overlap,  the  reduction  of 
compressive  stress  at  the  center  of  the  spot  is  negated  by  the  subsequent  shot. 
Another  interesting  result  is  that  at  all  the  locations  the  objective  is  achieved 
using  only  one  shot. 

5.4  Discussion 

The  optimization  of  structural  components  subjected  to  high-energy  impulses 
can  involve  nonlinear  elastic-plastic  behavior,  time-consuming  finite  element 
simulations,  and  mixed  optimization  variables.  To  solve  such  problems,  this 
research  develops  the  progressive  multifidelity  optimization  strategy.  The  strat¬ 
egy  begins  with  parametric  investigations  using  lower  fidelity  models  to  deter¬ 
mine  sensitive  parameters  and  to  localize  the  design  space.  In  each  step  of  the 
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strategy,  the  design  space  is  reduced  and  parameters  are  eliminated,  with  the 
goal  of  reaching  an  optimal  solution. 

The  demonstration  shows  that  pressure  pulse  magnitude,  spot  shape  and  size, 
and  layout  are  critical  parameters.  The  first  step  determines  the  shape  and  du¬ 
ration  of  the  pressure  pulse  that  is  suitable  for  the  formulated  problems.  In  the 
second  step,  a  bigger  spot  size  is  determined  better  for  the  given  formulation. 
Based  on  the  results  from  the  first  two  steps,  the  problem  is  formulated  for 
the  third  step.  Finally,  the  third  step  determines  a  five-shot  LP  layout  that  can 
meet  the  requirements  set  forth  in  the  optimization  formulation.  Compared  to 
a  traditional  approach,  the  combination  of  three  steps  and  the  use  of  surrogate 
models  substantially  reduces  the  computational  cost. 

Overall,  a  progressive  multifidelity  optimization  strategy  is  developed  in  this 
section.  A  surrogate  model  is  proposed  that  can  be  used  for  impulse-type 
simulations.  The  proposed  surrogate  model  accommodates  the  presence  of 
repeated  impulse  loads.  To  effectively  solve  the  optimization  problem,  mod¬ 
ifications  are  introduced  into  PSO:  LHS-based  particle  generation  and  a  dis¬ 
crete  variable  handling  technique.  Computational  time  is  further  reduced  by 
utilizing  the  parallel  processing  advantage  of  a  population-based  optimization 
technique.  The  strategy  is  validated  and  successfully  employed  to  determine 
the  optimal  LP  parameters.  This  design  strategy  and  the  proposed  surrogate 
model  can  be  applied  to  other  impulse-type  problems  with  mixed  variables  of 
optimization. 
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6.0  SUMMARY  AND  FUTURE  DIRECTIONS 


This  section  presents  the  contributions  and  summary  of  the  research  work  per¬ 
formed  in  this  report.  The  tasks  performed  in  this  report  involve  the  devel¬ 
opment  of  simulation  capabilities,  surrogate  models  for  a  flat  surface,  opti¬ 
mization  of  a  laser  peening  process  at  one-location,  and  optimization  of  a  laser 
peening  process  at  multiple-locations. 

6.1  Contributions 

6.1.1  EE  Simulation 

Modeling  and  Parameter  Design  of  a  Laser  Peening  Process.  [68, 113] 

A  framework  was  developed  for  parametric  simulations  that  can  consider  many 
of  the  LP  variables.  The  framework  incorporates  a  3D  model  drawn  from  the 
literature.  The  3D  model  was  used  to  investigate  the  LP  variables  such  as  tem¬ 
poral  variation  of  pressure,  pressure  pulse  magnitude,  and  spot  shape  and  size. 
Minor  modiflcations  to  the  symmetric  3D  model  allowed  us  to  investigate  com¬ 
ponent  thickness,  a  two-shot  sequence  at  multiple  locations,  and  a  seven-shot 
sequence  at  multiple  locations.  These  investigations  led  to  following  flndings: 

(i)  Spot  shape  can  be  employed  to  the  decrease  maximum  tensile  stress  with¬ 
out  signiflcantly  compromising  the  maximum  compressive  stress. 

(ii)  An  increase  in  pressure  magnitude  increases  compressive  depth,  volume, 
and  maximum  stress;  however,  compared  to  the  compressive  stress,  the 
increase  in  maximum  tensile  stress  is  greater. 

(hi)  The  depth  of  compressive  stress  tends  to  increase  with  thickness  and  sat¬ 
urates  after  5  mm.  The  maximum  tensile  stress  flrst  increases  and  then 
decreases  with  thickness. 

(iv)  Interaction  of  the  induced  residual  stress  proflle  of  two  sequential  shots 
is  dependent  upon  the  distance  between  the  two  shots.  The  interaction 
exists  even  when  there  is  no  overlap. 

(v)  In  the  considered  parameters  and  the  geometric  model  (symmetric  flat 
cylinder),  there  is  no  signiflcant  change  in  residual  stress  proflle  due  to  the 
sequence  of  shots;  however,  layout  and  number  of  shots  show  signiflcant 
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sensitivities. 


6.1.2  Optimization 
LP  Optimization:  One  location 

Modified  PSO  for  Multimodal  Mixed-variable  Optimization.  [114] 


A  modified  PSO  method  is  proposed  to  solve  a  multimodal  mixed-variable 
optimization  problem.  The  modified  PSO  employs  a  surrogate  model,  a  nich¬ 
ing  strategy,  and  a  mixed- variable  handling  technique:  the  surrogate  model 
is  employed  to  manage  computational  cost,  the  niching  strategy  is  used  to 
find  and  preserve  multiple  optima,  and  the  mixed- variable  technique  is  used 
to  manage  continuous  and  integer  variables.  The  method  is  demonstrated  on 
multiple  test  problems  and  on  a  one-location  laser  peening  optimization.  The 
solution  of  the  laser  peening  problem  provides  the  designer  with  multiple  op¬ 
tions.  The  designer  can  select  the  most  suitable  option,  depending  upon  the 
component  requirements. 

(i)  Multiple  solutions  exist  to  achieve  the  target  residual  stress  profile. 

(ii)  Modified  PSO  handles  mixed  variables  effectively  and  finds  multiple  op¬ 
timal  solutions. 

(iii)  Multiple  shots  at  lower  pressure  at  the  same  location  create  a  higher  depth 
of  compressive  stresses  compared  to  a  single  shot  of  higher  pressure. 

(iv)  One  shot  at  higher  pressure  creates  a  higher  magnitude  the  maximum 
compressive  stress  compared  to  multiple  shots  at  lower  pressure. 

(v)  An  elliptical  spot  shape,  compared  to  a  circular  spot  shape,  induces  a 
lower  magnitude  of  tensile  stress  for  similar  compressive  stress  magni¬ 
tude. 

Successful  demonstration  on  multiple  problems  show  that  the  proposed  method 
can  be  employed  to  solve  multimodal  mixed- variable  problems. 
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LP  Optimization:  Multiple  Location 


Mixed-variable  Optimization  Strategy  Employing  Multi-fidelity  Simulation  and 
Surrogate  Models.  [115,116] 

A  progressive,  simulation-based,  mixed- variable  optimization  strategy  is  de¬ 
veloped.  This  strategy  involves  three  steps  and  progressively  employs  axi- 
symmetric,  symmetric  3Z),  and  parametric  models.  A  surrogate  model  is  pro¬ 
posed  that  can  be  used  for  impulse-type  simulations.  The  proposed  surro¬ 
gate  models  accommodate  the  presence  of  repeated  impulse  loads.  To  effec¬ 
tively  solve  the  optimization  problem,  two  modifications  are  introduced  into 
PSO:  LHS-based  particle  generation  and  a  discrete  variable  handling  tech¬ 
nique.  Computational  time  is  further  reduced  by  utilizing  the  parallel  process¬ 
ing  advantage  of  a  population-based  optimization  technique.  The  strategy  is 
validated  and  successfully  employed  to  determine  the  optimal  LP  parameters. 
This  design  strategy  and  the  proposed  surrogate  model  can  be  applied  to  other 
impulse-type  problems  with  mixed  optimization  variables. 

(i)  A  progressive  three-step  optimization  strategy  was  developed  that  can  be 
applied  to  nonlinear  impulse  type  processes  such  as  laser  peening,  shot 
peening,  and  bullet  impacts  on  structures. 

(ii)  A  surrogate  model  is  proposed  to  alleviate  the  prohibitive  computational 
cost  of  an  iterative  optimization  process. 

(iii)  A  parallel  processing  technique  (MATLAB-i-PuTTy)  is  developed  for  ef¬ 
fective  population-based  optimization. 


6.2  Research  Summary 

Overall,  a  simulation  capability  is  developed  that  includes  a  symmetric  3D 
model,  a  parametric  plate,  and  rectangular  bar  coupon  models.  These  models 
were  employed  to  investigate  the  individual  LP  parameters,  optimize  single 
and  multiple  locations  peening.  The  investigations  and  optimization  assist 
an  engineer  to  achieve  desired  residual  stress  profile.  A  multimodal  mixed- 
variable  optimization  method  was  developed  to  solve  the  one-location  laser 
peening  optimization  problem.  The  method  modifies  particle  swarm  opti¬ 
mization  to  incorporate  multimodality  and  mixed-variable  handling  proper- 
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ties.  A  progressive  multifidelity  optimization  strategy  was  developed  to  solve 
the  multiple-location  laser  peening  optimization  problem.  The  PMOS  com¬ 
bines  modified  PSO,  LHS,  mixed- variable  strategy,  and  multifidelity  surrogate 
models  to  solve  a  laser  peening  optimization  problem. 

6.3  Future  Direction 

The  presented  research  work  can  be  continued  in  many  directions.  The  four 
major  objectives  can  be  reduction  in  the  simulation  time,  parametric  interac¬ 
tions  effects,  process  optimization  for  a  practical  geometry,  and  uncertainty 
quantification  of  the  LP  process.  The  parametric  investigations,  optimization, 
and  uncertainty  quantification  can  be  simplified  if  the  simulation  time  can  be 
reduced  by  further  research.  Empirical  formulations,  surrogate  models,  and 
spectral  finite  element  method  [1171  can  be  investigated  for  this  purpose.  The 
parametric  investigations  performed  in  this  research  are  limited  to  one  variable 
at  a  time,  ignoring  the  interaction  effects.  It  will  be  interesting  to  determine 
various  interaction  effects  of  the  parameters. 

Determining  optimized  process  parameters  for  a  practical  3Z)  problem  can  be 
very  useful  for  the  industrial  application.  This  is  possible  by  aggressively  em¬ 
ploying  surrogate  model.  Surrogate  models  developed  in  this  research  can  be 
advanced  to  use  in  a  practical  problem  such  as  aircraft  lug.  As  discussed  in  the 
report,  LP  is  a  mixed-variable  process.  The  development  of  a  mixed-variable 
uncertainty  quantification  method  can  assist  in  uncertainty  quantification  of 
the  LP  process.  Additional  potential  areas  include  experimental  validation  of 
multiple-location  peening,  optimization,  and  fatigue  life  optimization  strate¬ 
gies  [1 181.  This  experimental  work  will  not  only  provide  additional  validation 
of  the  simulation  work  but  also  facilitate  the  development  of  an  optimization 
strategy  for  practical  surface  enhancement  and  fatigue  problems. 
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Ad 

£() 

^eq 

<3hEL 

^surf 

^yd 

A 

ACO 

B 

C 

Cd 

CD 

CV 

d 

d* 

Dp 

DOE 

EV 

EEA 

GA 

GBP 

HEL 

/o 

Le 

LASER 

LHS 

LP 

LPB 

m 

mbesti 

MC 

MGPSO 

MMVO 

MNPSO 


Change  in  Depth 
Normalized  Plastic  Strain  Rate 
Effective  Stress 
Hugoniot  Elastic  Limit 
Surface  Residual  Stress 
Dynamic  Yield  Stress 

Material  Constant  for  the  Johnson-Cook  Material  Model 
Ant  Colony  Optimization 

Material  Constant  for  the  Johnson-Cook  Material  Model 
Material  Constant  for  the  Johnson-Cook  Material  Model 
Current  Effective  Dilatational  Wave  Speed  of  the  Material 
Depth  of  Compressive  Stress 
Total  Compressive  Volume 
Depth 

Inclusive  Distance  Surrounding  mbesti 

Plastically  Affected  Depth 

Design  of  Experiments 

Evolutionary  Algorithms 

Einite  Element  Analysis 

Genetic  Algorithms 

Glass  Bead  Peening 

Hugoniot  Elastic  Limit 

Incident  Laser  Power  Density 

Minimum  Element  Size 

Light  Amplification  by  Simulated  Emission  of  Radiation 

Latin  Hypercube  Sampling 

Laser  Peening 

Low  Plasticity  Burnishing 

Material  Constant  for  the  Johnson-Cook  Material  Model 

Global  or  Local  Best  Eound  so  Ear 

Maximum  Compressive  Stress 

Multigrouped  Particle  Swarm  Optimization 

Multimodal  Mixed- Variable  Optimization 

Mixed- variable  Niche  Particle  Swarm  Optimization 
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MT 

Maximum  Tensile  Stress 

MVO 

Mixed  Variable  Optimization 

OSC 

Ohio  Super-computing  Center 

P 

Pressure  Pulse  Magnitude 

P 

Pressure 

PCX 

Parallel  Computing  Toolbox 

PMOS 

Progressive  Multifidelity  Optimization  Strategy 

PSO 

Particle  Swarm  Optimization 

q 

Specified  Number  of  Particles  Allowed  to  a  Single  Optima 

r 

Shot  Radius 

RSM 

Response  Surface  Method 

Spot  Shape 

SA 

Simulated  Annealing 

SP 

Shot  Peening 

t 

Pressure  Pulse  Midspan 

Tjnelt 

Melting  Temperature 

TANA 

Two-point  Adaptive  Nonlinear  Approximation 

TANA2 

Improved  Two-point  Adaptive  Nonlinear  Approximation 

V. 

Velocity  of  Elastic  Wave 

Velocity  of  Plastic  Wave 

VISAR 

Velocity  Interferometer  System  for  any  Reflector 

WP 

Water  Peening 
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